/*
 * transient_search
 * Jamie Stevens
 *
 * This program will go through a bunch of RPFITS files to assess the data they contain.
 *
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "RPFITS.h"
#include <sys/types.h>
#include <sys/resource.h>
#include <dirent.h>
#include <fnmatch.h>
#include <time.h>
#include <sys/time.h>
#include "cpgplot.h"
#include <complex.h>
#include <fftw3.h>
#include <argp.h>
#include <gsl/gsl_fit.h>
#include <signal.h>
#include <libxml/encoding.h>
#include <libxml/xmlwriter.h>
#include <libxml/xmlreader.h>

/* the default size for strings */
#define BUFSIZE 256

#define YES 1
#define NO  0

/* XML encoding type */
#define MY_ENCODING "ISO-8859-1"

/* RPFITS commands numerical definitions to make the code more readable */
#define JSTAT_OPENFILE -3
#define JSTAT_OPENFILE_READHEADER -2
#define JSTAT_READNEXTHEADER -1
#define JSTAT_READDATA 0
#define JSTAT_CLOSEFILE 1
#define JSTAT_SKIPTOEND 2

/* RPFITS return codes definitions to make the code more readable */
#define JSTAT_UNSUCCESSFUL -1
#define JSTAT_SUCCESSFUL 0
#define JSTAT_HEADERNOTDATA 1
#define JSTAT_ENDOFSCAN 2
#define JSTAT_ENDOFFILE 3
#define JSTAT_FGTABLE 4
#define JSTAT_ILLEGALDATA 5

/* matching criteria are set here and checked programmatically 
   making it easier to tweak what we're looking for */
/* the number of instruments that can be considered */
#define CONSIDER_N_INSTRUMENTS 1
/* the list of instruments that can be considered */
char consider_instrument[CONSIDER_N_INSTRUMENTS][BUFSIZE]={"ATCA"};

/* the minimum number of IFs required */
#define MIN_N_IFS 2
/* uncomment this next line if this is a minimum number instead of an exact
   number */
/* #define MORE_IFS */

/* uncomment this next line if you don't care about simultaneity */
/* #define NOT_SIMULTANEOUS */

/* the maximum number of frequency channels */
#define MAX_N_CHANNELS 32
/* uncomment this next line if this is a maximum number instead of an exact
   number */
/* #define LESS_CHANNELS */

/* the maximum frequency to consider, in Hz */
#define MAX_IF_FREQUENCY 1600E6

/* the minimum bandwidth, in Hz */
#define MIN_BANDWIDTH 64E6

/* don't count visibilities within this number of from the scan boundaries */
#define SCAN_BOUNDARY_LIMIT 5
/* uncomment this line if you want to keep data at the end of scans */
#define KEEP_ENDOFSCAN

/* the different types of plot that can be produced */
#define PLOT_INDIVIDUAL_SPECTRA 2
#define PLOT_WATERFALL 4
#define PLOT_ALL 8
#define PLOT_SEPARATE 16
#define PLOT_COMBINED_POLS 32
#define PLOT_PEAKS 64
#define PLOT_NONE 128

/* the different colour tables we can use in PGPLOT */
#define PGPLOT_CT_GREYSCALE 1
#define PGPLOT_CT_RAINBOW 2
#define PGPLOT_CT_HEAT 3
#define PGPLOT_CT_WEIRDIRAF 4
#define PGPLOT_CT_AIPS 5

/* the different Stokes parameters we can make */
#define STOKES_I -1
#define STOKES_Q -2
#define STOKES_U -3
#define STOKES_V -4

/* some PGPLOT constants that I should really put into a global header
   file for future use */
#define PGSCI_BACKGROUND 0
#define PGSCI_INVERSE 1
#define PGSCI_RED 2
#define PGSCI_GREEN 3
#define PGSCI_BLUE 4
#define PGSCI_CYAN 5
#define PGSCI_MAGENTA 6
#define PGSCI_YELLOW 7
#define PGSCI_ORANGE 8
#define PGSCI_GREENYELLOW 9
#define PGSCI_GREENCYAN 10
#define PGSCI_BLUECYAN 11
#define PGSCI_BLUEMAGENTA 12
#define PGSCI_REDMAGENTA 13
#define PGSCI_DARKGREY 14
#define PGSCI_LIGHTGREY 15

#define PGSLS_FULL 1
#define PGSLS_DASHED 2
#define PGSLS_DOTDASH 3
#define PGSLS_DOTTED 4
#define PGSLS_DASHDOTS 5

#define PGLEN_NORMALIZED_DEVICE 0
#define PGLEN_INCHES 1
#define PGLEN_MM 2
#define PGLEN_ABSOLUTE_DEVICE 3
#define PGLEN_WORLD 4
#define PGLEN_VIEWPORT_FRACTION 5

#define PGQINF_ITEM_VERSION "VERSION"
#define PGQINF_ITEM_STATE "STATE"
#define PGQINF_ITEM_USER "USER"
#define PGQINF_ITEM_DATETIME "NOW"
#define PGQINF_ITEM_DEVICE "DEVICE"
#define PGQINF_ITEM_TYPE "TYPE"
#define PGQINF_ITEM_DEVICETYPE "DEV/TYPE"
#define PGQINF_ITEM_ISHARDCOPY "HARDCOPY"
#define PGQINF_ITEM_ISTERMINAL "TERMINAL"
#define PGQINF_ITEM_HASCURSOR "CURSOR"
#define PGQINF_ITEM_HASSCROLL "SCROLL"
#define PGQINF_VALUE_UNSUCCESFUL "?"
#define PGQINF_VALUE_YES "YES"
#define PGQINF_VALUE_NO "NO"
#define PGQINF_VALUE_OPEN "OPEN"
#define PGQINF_VALUE_CLOSED "CLOSED"

#define PEAK_VERIFIED_NO -1
#define PEAK_VERIFIED_FIRST_STAGE -2
#define PEAK_VERIFIED_YES -3

typedef struct spectral_info_ {
  float average_amplitude;
  float rms_noise;
  int npeaks_found;
  int *peak_channels;
  int *peak_verified;
  int *peak_number;
  int *peak_present_both_pols;
  float estimated_average_amplitude;
  float estimated_average_amplitude_error;
  float estimated_rms_noise;
  double fit_c0;
  double fit_c1;
} spectral_info;

/* the structure to hold all the data in */
typedef struct correlation_data_ {
  float *vis;
  float *wgt;
  int baseline;
  int islongest_baseline;
  float ut;
  time_t observation_time;
  float u;
  float v;
  float w;
  int flag;
  int bin;
  int if_no;
  int sourceno;
  char sourcename[BUFSIZE];
  char **stokes_name;
  int if_nfreq;
  int if_nstok;
  float if_bandwidth;
  int scannum;
  float *pols_combined;
  float **individual_pols;
/*   int polarisation_non_detection; */
  char combined_pol_name[BUFSIZE];
  spectral_info spectrum_info;
  struct correlation_data_ *next;
  struct correlation_data_ *prev;
  struct correlation_data_ *separate_next;
  struct correlation_data_ *samebaseline_next;
  struct correlation_data_ *samebaseline_prev;
} correlation_data;

typedef struct peak_baseline_ {
  int if_number;
  int baseline_ant1;
  int baseline_ant2;
  float baseline_length;
  int lag_channel;
  float baseline_lag;
  float baseline_lag_err;
  float amplitude;
  float spectrum_average_amplitude;
  float spectrum_rms_noise;
  float significance;
} peak_baseline;

typedef struct transient_peak_ {
  int peak_number;
  char source_name[BUFSIZE];
  time_t observation_time;
  int n_lags;
  peak_baseline **lag_peaks;
  int available_baselines;
  int peak_baselines;
  float peak_baselines_fraction;
  int transient_candidate;
} transient_peak;

typedef struct archive_file_ {
  char path[BUFSIZE];
  char filename[BUFSIZE];
  int usable;
  int n_peaks;
  transient_peak **transient_peaks;
  char plot_name[BUFSIZE];
  struct archive_file_ *next;
} archive_file;

int opened_pgplot_device_id;
char opened_pgplot_device_name[BUFSIZE];

/* routine definitions */
void strip_path(char *full_path,char *filename);
void catch_interrupt(int sig);
void stop_program(int returncode);
void strip_trailing_spaces(char *string);
int choose_CA_files(const struct dirent *ep);
void baseline_to_antennae(int baseline_number,int *antenna1,int *antenna2);
int want_this_file(int *longest_baseline);
int read_data(correlation_data **all_correlations,int longest_baseline,int *keep_reading);
int output_medium(char *PGPLOT_device);
int plot_data(correlation_data **all_correlations,int plot_type);
void which_panel(int baseline,int if_no,int *xpanel,int *ypanel);
int add_correlation(correlation_data **correlation_list,float **vis,float **wgt,int baseline,
		    float ut,float u,float v,float w,int flag,int bin,int if_no,int sourceno,
		    int if_nfreq,int if_nstok,int longest_baseline,int scannum);
int clear_correlations(correlation_data **all_correlations);
int remove_correlation(correlation_data **correlation_to_remove);
int separate_lists(correlation_data **all_correlations,
		   correlation_data ***separate_correlations);
void palett(int type,float contra,float bright);
void setvp(void);
int change_to_lagspace(correlation_data **all_correlations,
		       correlation_data **all_transforms,
		       int cycle_type);
void copy_correlation_info(correlation_data **original,
			   correlation_data **destination);
int combine_polarisations(correlation_data **all_correlations,int cycle_type,int stokes_type);
int peak_finder(correlation_data **all_transforms,correlation_data ***peak_lists);
int cull_peaks(correlation_data **all_transforms);
int check_for_overlap(float value1,float value1_err,float value2,float value2_err);
int time_sort_data(correlation_data **data);
float standard_deviation(float **values,int nelements,float *average);
float meanvalue(float **values,int nelements);
float standard_deviation_exclude(float **values,int nelements,float *average,
				 int n_exclusions,int **exclusions);
archive_file *file_info(char *filename);
void build_file_database(void);
xmlChar * ConvertInput(const char *in, const char *encoding);
void lag_calculator(int channel_number,float bandwidth,float *lag,float *lag_err);
int channel_calculator(float bandwidth,float lag_value,float lag_err);
void lag_scaler(float bandwidth,int long_baseline_channel,
		float long_baseline_u,float long_baseline_v,float long_baseline_w,
		float this_baseline_u,float this_baseline_v,float this_baseline_w,
		float *this_baseline_lag,float *this_baseline_lag_err,
		int *this_baseline_channel);
void additional_peak_checks(correlation_data *transform,int peak_number);
void shift_data(float *data,int ndata,int nshift);
void findpeaks(float *data,int ndata,float slopethreshold,float ampthreshold,
	       float smoothwidth,int peakgroup,int *npeaks,int **peakchannel,
	       float **peakheight,float **peakwidth);
void fastsmooth(float *data,int ndata,int smoothwidth,float **smoothed_data);
void deriv(float *data,int ndata,float **deriv_data);
void val2ind(float *data,int ndata,float val,int **indices,float **values,
	     int *nindices);

/* now the argp command-line argument parsing stuff */
static struct argp_option options[] = {
  {"debug",'d',0,0,"print debugging information"},
  {"pgplot-device",'P',"PGPLOT_DEVICE",0,"the PGPLOT device to use for graphics output"},
  {"pgplot-ask",'A',0,0,"make PGPLOT ask before erasing page"},
  {"plot-type",'p',"PLOT_TYPE",0,"the type of plot to make (waterfall/spectra/peaks/none)"},
  {"nearby-spectra",'n',"NEARBY_SPECTRA",0,"the number of nearby spectra to use when searching for peaks"},
  {"sigma-factor",'s',"SIGMA_FACTOR",0,"how many standard deviations must the peak be to be detected"},
  {"report-file",'r',"REPORT_FILE",0,"make a full report to this filename"},
  {"files-file",'f',"FILES_FILE",0,"specify the list of files to look at"},
  {"min-match-fraction",'m',"MIN_MATCH_FRACTION",0,"the minimum fraction of baselines that must have a peak"},
  {"existing-xml-report",'x',"EXISTING_XML_REPORT",0,"name of an existing XML report"},
  { 0 }
};

struct arguments {
  int DEBUGGING;
  char PGPLOT_device[BUFSIZE];
  int PGPLOT_ask;
  int plot_type;
  int number_of_nearby_spectra;
  float minimum_peak_sigma_factor;
  char report_filename[BUFSIZE];
  char xml_report_filename[BUFSIZE];
  int make_report;
  int make_xml_report;
  FILE *report_fh;
  char files_filename[BUFSIZE];
  float minimum_matching_fraction;
  xmlTextWriterPtr xml_report_writer;
  int use_existing_xml_report;
  char existing_xml_report[BUFSIZE];
  archive_file *existing_xml_report_files;
};

struct arguments arguments;

char report_error[BUFSIZE];

/* some macros for reporting status, errors and debugging messages */
#define REPORT(...) if (arguments.make_report==1) fprintf(arguments.report_fh, __VA_ARGS__);
#define DEBUG(...) if (arguments.DEBUGGING==1){ printf("++ "); printf(__VA_ARGS__); }
#define STATUS(...) printf("== "); printf(__VA_ARGS__);
#define ERROR(...) fprintf(stderr,"!! "); fprintf(stderr, __VA_ARGS__); snprintf(report_error,BUFSIZE, __VA_ARGS__);
#define XMLREPORT(routine,...) if (arguments.make_xml_report==1){ int rc; rc=routine(__VA_ARGS__); if (rc<0) { ERROR("XML writing error at %s:%d\n",__FILE__,__LINE__); arguments.make_xml_report=0; } }
#define XMLCONVERT(string,xmlstring,...) if (xmlstring!=NULL){ xmlFree(xmlstring); } snprintf(string,BUFSIZE,__VA_ARGS__); xmlstring=ConvertInput(string,MY_ENCODING);

/* make a safe malloc, realloc, calloc & free */
#define MALLOC(a,b) if ((a=malloc((b)))==NULL) { ERROR("malloc failed at %s:%d!\n",__FILE__,__LINE__); stop_program(-1); }
#define REALLOC(a,b) if ((a=realloc(a,(b)))==NULL) { ERROR("realloc failed at %s:%d!\n",__FILE__,__LINE__); stop_program(-1); }
#define CALLOC(a,b,c) if ((a=calloc(b,c))==NULL) { ERROR("calloc failed at %s:%d!\n",__FILE__,__LINE__); stop_program(-1); }
#define FREE(a) { free(a); a=NULL; }

static error_t parse_opt(int key,char *arg,struct argp_state *state){
  struct arguments *arguments = state->input;

  switch (key){
  case 'd': /* debugging output */
    arguments->DEBUGGING=1;
    STATUS("debugging output enabled\n");
    break;
  case 'P': /* PGPLOT device */
    strcpy(arguments->PGPLOT_device,arg);
    STATUS("will output plot to the PGPLOT device %s\n",arguments->PGPLOT_device);
    break;
  case 'A': /* make PGPLOT ask before erasing page */
    arguments->PGPLOT_ask=1;
    STATUS("PGPLOT will ask before moving to a new page\n");
    break;
  case 'p': /* type of plot to create */
    if (strcmp(arg,"waterfall")==0)
      arguments->plot_type=PLOT_WATERFALL;
    else if (strcmp(arg,"spectra")==0)
      arguments->plot_type=PLOT_INDIVIDUAL_SPECTRA;
    else if (strcmp(arg,"peaks")==0)
      arguments->plot_type=PLOT_PEAKS;
    else if (strcmp(arg,"none")==0)
      arguments->plot_type=PLOT_NONE;
    else {
      /* unknown plot type */
      snprintf(report_error,BUFSIZE,"Do not recognise plot type %s\n",arg);
      ERROR(report_error);
      stop_program(-1);
    }
    STATUS("will create a %s plot\n",arg);
    break;
  case 'n': /* the number of nearby spectra */
    arguments->number_of_nearby_spectra=atoi(arg);
    STATUS("using +/- %d nearby spectra\n",arguments->number_of_nearby_spectra);
    break;
  case 's': /* the "sigma factor" */
    arguments->minimum_peak_sigma_factor=atof(arg);
    STATUS("setting the minimum peak sigma factor to %.1f\n",arguments->minimum_peak_sigma_factor);
    break;
  case 'r': /* make a report file */
    arguments->make_report=1;
    arguments->make_xml_report=1;
    strncpy(arguments->report_filename,arg,BUFSIZE);
    strncpy(arguments->xml_report_filename,arg,BUFSIZE);
    strcat(arguments->xml_report_filename,".xml");
    STATUS("will output full report to file %s\n",arguments->report_filename);
    STATUS("will output XML report to file %s\n",arguments->xml_report_filename);
    break;
  case 'f': /* a file containing the list of RPFITS files to read */
    strncpy(arguments->files_filename,arg,BUFSIZE);
    STATUS("will get files to read from file %s\n",arguments->files_filename);
    break;
  case 'm': /* the minimum fraction of baselines with a corresponding peak
	       needed to consider the peak a transient */
    arguments->minimum_matching_fraction=atof(arg);
    STATUS("setting minimum match fraction to %.3f\n",arguments->minimum_matching_fraction);
    break;
  case 'x': /* an old XML report file that we will consult to prevent
	       unnecessary effort */
    arguments->use_existing_xml_report=1;
    strncpy(arguments->existing_xml_report,arg,BUFSIZE);
    STATUS("will use existing XML report file %s\n",arguments->existing_xml_report);
    break;
  }
  
  return(0);
}

static struct argp argp = { options, parse_opt, 0, 0 };

void strip_path(char *full_path,char *filename){
  /* this routine takes a full path name of a file and returns just the
     filename without the path information */
  int i;

  /* is there a slash in this filename? */
  if (strstr(full_path,"/")==NULL){
    /* no slash, so we just copy the whole string */
    strcpy(filename,full_path);
  } else {
    /* there is a slash */
    for (i=strlen(full_path)-1;i>=0;i--){
      if (full_path[i]=='/'){
	/* this is the last slash */
	strcpy(filename,full_path+i+1);
	break;
      }
    }
  }
}

void catch_interrupt(int sig){
  /* this routine should catch the Ctrl-C signal and call the stop_program routine
     so that the program cleans up properly */
  STATUS("\n");
  if (sig==SIGINT){
    STATUS("received interrupt signal\n");
    stop_program(-2);
  }
}

void stop_program(int returncode){
  /* this routine is called whenever the program needs to end, so it needs to do 
     some cleanup and such to make things neat */

  /* if we're writing a report we should write out why we're stopping and then 
     close the file */
  REPORT(report_error);
  if (arguments.make_report==1)
    fclose(arguments.report_fh);

  /* if we're writing an XML report we should end the document */
  XMLREPORT(xmlTextWriterEndDocument,arguments.xml_report_writer);
  if (arguments.make_xml_report==1)
    xmlFreeTextWriter(arguments.xml_report_writer);

  /* close any open PGPLOT devices */
  if (strlen(opened_pgplot_device_name)!=0){
    cpgclos();
    opened_pgplot_device_name[0]='\0';
  }

  STATUS("Program ends.\n");
  exit(returncode);
}

void strip_trailing_spaces(char *string){
  /* this routine takes a string that may or may not have trailing spaces before the null 
     terminator, and removes those trailing spaces if they exist */
  int i;
  
  for (i=strlen(string)-1;i>=0;i--){
    if (string[i]==' '){
      string[i]='\0'; /* replace the space with the null terminator */
    } else {
      /* not a space so we stop here */
      break;
    }
  }
}

/* some macros for making benchmark measurements of program execution time */
#define START_TIMER gettimeofday(&start_benchmark,NULL)
#define STOP_TIMER gettimeofday(&end_benchmark,NULL)
#define TIMER_BENCHMARK benchmark_time=(float)(end_benchmark.tv_sec-start_benchmark.tv_sec)+(float)(end_benchmark.tv_usec-start_benchmark.tv_usec)/1E6;
#define PRINT_BENCHMARK(a) if (benchmark_time>=10.0) { snprintf(a,BUFSIZE,"%.5f secs",benchmark_time); } else { snprintf(a,BUFSIZE,"%.2f msecs",benchmark_time*1E3); }
#define START_COUNTER getrusage(RUSAGE_SELF,&start_rusage);
#define STOP_COUNTER getrusage(RUSAGE_SELF,&end_rusage);
#define COUNTER_BENCHMARK blocks_read=end_rusage.ru_inblock-start_rusage.ru_inblock;

int main(int argc,char *argv[]){
  
  int i,j,rpfits_result,wtf_result,n_lists;
  int this_jstat,keep_reading,longest_baseline;
  float *vis=NULL,*wgt=NULL,ut,u,v,w,added_time=0,benchmark_time;
  float total_file_time=0,total_total_time=0;
  int baseline,flag,bin,if_no,sourceno,usable_data;
  int ngood=0,header_number,correlations_added=0,blocks_read;
  FILE *good_files=NULL;
  char fline[BUFSIZE],**goodfiles=NULL,files_file[BUFSIZE];
  char printed_benchmark[BUFSIZE],tmp[BUFSIZE];
  correlation_data *all_correlations=NULL,**separate_correlations=NULL;
  correlation_data *all_transforms=NULL,**found_peaks=NULL;
  time_t time_now;
  struct tm now_time;
  struct timeval start_benchmark,end_benchmark;
  struct rusage start_rusage,end_rusage;
  xmlChar *xmltmp=NULL;
  archive_file *existing_file=NULL;

  /* initialise the arguments and some other values */
  arguments.DEBUGGING=0;
  arguments.PGPLOT_ask=0;
  arguments.plot_type=PLOT_INDIVIDUAL_SPECTRA;
  arguments.number_of_nearby_spectra=3;
  arguments.minimum_peak_sigma_factor=3.0;
  strcpy(arguments.PGPLOT_device,"/xs");
  arguments.report_filename[0]='\0';
  arguments.xml_report_filename[0]='\0';
  arguments.make_report=0;
  arguments.make_xml_report=0;
  arguments.report_fh=NULL;
  report_error[0]='\0';
  arguments.files_filename[0]='\0';
  opened_pgplot_device_name[0]='\0';
  opened_pgplot_device_id=-1;
  arguments.minimum_matching_fraction=0.90; /* 90% of baselines need the same peak */

  STATUS("TRANSIENT SEARCH PROGRAM\n");
  time(&time_now);
  gmtime_r(&time_now,&now_time);
  STATUS("Program started at %02d:%02d:%02d %02d/%02d/%04d GMT\n",
	 now_time.tm_hour,now_time.tm_min,now_time.tm_sec,
	 now_time.tm_mday,now_time.tm_mon+1,now_time.tm_year+1900);
  STATUS("\n");

  /* now parse the command line arguments */
  argp_parse(&argp,argc,argv,0,0,&arguments);
  STATUS("\n");

  /* now setup the signal handler to catch interrupts */
  signal(SIGINT,catch_interrupt);
  /* and to ignore any hangups */
  signal(SIGHUP,SIG_IGN);

  /* are we making a report? */
  if (arguments.make_report==1){
    if ((arguments.report_fh=fopen(arguments.report_filename,"w"))==NULL){
      ERROR("unable to open report file %s for writing! will not create a report!\n",
	     arguments.report_filename);
      arguments.make_report=0;
    }
  }
  /* we also try to make an XML report as it should be easier to quickly get data
     from later */
  if (arguments.make_xml_report==1){
    arguments.xml_report_writer=xmlNewTextWriterFilename(arguments.xml_report_filename,0);
    if (arguments.xml_report_writer==NULL){
      ERROR("unable to open XML report file %s for writing! will not create an XML report!\n",
	    arguments.xml_report_filename);
      arguments.make_xml_report=0;
    }
  }
  /* are we going to consulting an exisiting XML file? */
  if (arguments.use_existing_xml_report==1){
    build_file_database();
  }

  /* the stuff up the top of the report file */
  REPORT("TRANSIENT SEARCH PROGRAM\n");
  REPORT("Program started at %02d:%02d:%02d %02d/%02d/%04d GMT\n",
	 now_time.tm_hour,now_time.tm_min,now_time.tm_sec,
	 now_time.tm_mday,now_time.tm_mon+1,now_time.tm_year+1900);
  REPORT("\n");

  REPORT("Will estimate average power and RMS noise using +/- %d nearby spectra\n",
	 arguments.number_of_nearby_spectra);
  REPORT("Will identify peaks that are %.1f RMS noise levels above the average level\n",
	 arguments.minimum_peak_sigma_factor);
  REPORT("%.1f %% of baselines must have a corresponding peak for transient detection\n",
	 arguments.minimum_matching_fraction*100);
  REPORT("\n");

  /* start the XML report with the XML default for the version, encoding and standalone
     declaration */
  XMLREPORT(xmlTextWriterStartDocument,arguments.xml_report_writer,NULL,MY_ENCODING,NULL);
  /* start the root element, call it ATCAARCHIVE */
  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "ATCAARCHIVE");
  /* start an element called DATETIME as a child of ATCAARCHIVE */
  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "DATETIME");
  /* now give it the attribute date and value */
  snprintf(tmp,BUFSIZE,"%02d/%02d/%04d",now_time.tm_mday,now_time.tm_mon+1,
	   now_time.tm_year+1900);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,BAD_CAST "date",
	    BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%02d:%02d:%02d",now_time.tm_hour,now_time.tm_min,now_time.tm_sec);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,BAD_CAST "time",
	    BAD_CAST tmp);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,BAD_CAST "timezone",
	    BAD_CAST "GMT");
  /* comments must come after attributes */
  XMLCONVERT(tmp,xmltmp,"the date and time at program start");
  XMLREPORT(xmlTextWriterWriteComment,arguments.xml_report_writer,xmltmp);
  /* close the element DATETIME */
  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
  /* now write out the parameters we're using */
  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "PARAMETERS");
  snprintf(tmp,BUFSIZE,"%d",arguments.number_of_nearby_spectra);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "nearby_spectra",BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%.1f",arguments.minimum_peak_sigma_factor);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "minimum_significance", BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%.3f",arguments.minimum_matching_fraction);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "baseline_match_fraction",BAD_CAST tmp);
  XMLCONVERT(tmp,xmltmp,"nearby_spectra is the number of nearby spectra use to estimate the average power and RMS noise of a target spectrum");
  XMLREPORT(xmlTextWriterWriteComment,arguments.xml_report_writer,xmltmp);
  XMLCONVERT(tmp,xmltmp,"minimum_significance is the minimum number of RMS noise levels above the average a peak needs to be to be considered");
  XMLREPORT(xmlTextWriterWriteComment,arguments.xml_report_writer,xmltmp);
  XMLCONVERT(tmp,xmltmp,"baseline_match_fraction is the fraction of baselines that must have a corresponding peak to be considered");
  XMLREPORT(xmlTextWriterWriteComment,arguments.xml_report_writer,xmltmp);
  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	     
  /* open a file which tells us which RPFITS files to look at */
  if (strlen(arguments.files_filename)==0){
#ifdef ATNF
    strcpy(files_file,"/u/ste616/transient_search/accepted_files");
#endif
#ifdef SEBBEN
    strcpy(files_file,"./accepted_files");
#endif
  } else {
    strncpy(files_file,arguments.files_filename,BUFSIZE);
  }
  STATUS("getting list of RPFITS files from file %s\n",files_file);
  if ((good_files=fopen(files_file,"r"))==NULL){
    ERROR("Can't open file %s for output of good files!\n",
	  files_file);
    stop_program(-1);
  }

  /* read the list of files that have good data */
  while((fgets(fline,BUFSIZE,good_files))!=NULL){
    (strstr(fline,"\n"))[0]='\0';
    if ((fline[0]!='-')&&(fline[0]!='+')){
      ngood++;
      REALLOC(goodfiles,ngood*sizeof(char*));
      MALLOC(goodfiles[ngood-1],BUFSIZE*sizeof(char*));
      strcpy(goodfiles[ngood-1],fline);
    }
  }
  fclose(good_files);

  total_total_time=0;
  for (i=0;i<ngood;i++){
    total_file_time=0;
    usable_data=0;
    STATUS("found RPFITS file %s\n",goodfiles[i]);
    /* check whether the existing XML file has data on this already */
    if (arguments.use_existing_xml_report==1){
      existing_file=file_info(goodfiles[i]);
      if (existing_file!=NULL){
	STATUS("have previously searched this file for transients, skipping\n");
	continue;
      }
    }
    strcpy(names_.file,goodfiles[i]); /* this tells the RPFITS routines which file to open */
    START_COUNTER;
    START_TIMER;
    this_jstat=JSTAT_OPENFILE;
    rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,&flag,&bin,&if_no,&sourceno);
    if (this_jstat==JSTAT_UNSUCCESSFUL){
      /* we can't open this file for some reason */
      ERROR("Unable to open RPFITS file %s\n",names_.file);
      REPORT(report_error);
      continue; /* go to the next file */
    }
    STATUS("reading RPFITS file\n");
    REPORT("Now working on RPFITS file %s\n",goodfiles[i]);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "FILE");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "filename",BAD_CAST goodfiles[i]);
    keep_reading=1;
    header_number=0;
    while(keep_reading){
      /* read the next header */
      this_jstat=JSTAT_READNEXTHEADER;
      rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,&flag,&bin,&if_no,&sourceno);
      DEBUG("rpfits_result = %d, baseline = %d\n",this_jstat,baseline);
      DEBUG("flag = %d, if_no = %d, sourceno = %d\n",flag,if_no,sourceno);
      if (this_jstat==JSTAT_SUCCESSFUL){
	header_number++;
	REPORT("Found header %d:\n",header_number);
	XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "HEADER");
	snprintf(tmp,BUFSIZE,"%d",header_number);
	XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		  BAD_CAST "number",BAD_CAST tmp);
	/* figure out if this data is useful */
	wtf_result=want_this_file(&longest_baseline);
	if (wtf_result==1){
	  /* yes it is useful */
	  usable_data=1;
	  DEBUG("we will use this data\n");
	  correlations_added=read_data(&all_correlations,longest_baseline,&keep_reading);
	  added_time=(float)correlations_added*(float)param_.intime;
	  REPORT("This %.1f seconds of data will be searched for transients.\n",
		 added_time);
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SUMMARY");
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "usable",BAD_CAST "yes");
	  snprintf(tmp,BUFSIZE,"%.1f",added_time);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "usable_time",BAD_CAST tmp);
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	  total_file_time+=added_time;
	  total_total_time+=added_time;
	} else {
	  DEBUG("we will not use this data\n");
	  REPORT("This data is unsuitable for our purposes.\n");
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SUMMARY");
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "usable",BAD_CAST "no");
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	}
	XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
      } else if (this_jstat==JSTAT_ENDOFFILE){
	keep_reading=0;
      }
    }
    /* close the file */
    this_jstat=JSTAT_CLOSEFILE;
    rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,&flag,&bin,&if_no,&sourceno);
    STOP_TIMER;
    STOP_COUNTER;
    TIMER_BENCHMARK;
    COUNTER_BENCHMARK;
    PRINT_BENCHMARK(printed_benchmark);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "file_read_time",BAD_CAST printed_benchmark);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    if (usable_data==1){
      STATUS("there is %.1f seconds of usable data in this file\n",total_file_time);
      XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SUMMARY");
      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		BAD_CAST "usable",BAD_CAST "yes");
      snprintf(tmp,BUFSIZE,"%.1f",total_file_time);
      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		BAD_CAST "total_usable_time",BAD_CAST tmp);
      XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    } else {
      STATUS("there is no usable data in this file\n");
      XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SUMMARY");
      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		BAD_CAST "usable",BAD_CAST "no");
      XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    }
    REPORT("Read %d blocks from RPFITS file %s in %s (%.2f blocks/sec).\n",
	   blocks_read,goodfiles[i],printed_benchmark,
	   ((float)blocks_read/benchmark_time));
    REPORT("Have so far read a total of %.1f seconds of usable data.\n",
	   total_total_time);

    /* take the Fourier transform of all the data */
    STATUS("computing Fourier transforms\n");
    START_TIMER;
    change_to_lagspace(&all_correlations,&all_transforms,PLOT_ALL);
    STOP_TIMER;
    TIMER_BENCHMARK;
    PRINT_BENCHMARK(printed_benchmark);
    REPORT("Fourier transforms took %s.\n",printed_benchmark);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "fourier_transform_time",BAD_CAST printed_benchmark);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);

    /* so now don't need all_correlations anymore */
    clear_correlations(&all_correlations);

    /* scalar add the polarisations */
    STATUS("combining polarisations\n");
    START_TIMER;
    combine_polarisations(&all_transforms,PLOT_ALL,STOKES_I);
    STOP_TIMER;
    TIMER_BENCHMARK;
    PRINT_BENCHMARK(printed_benchmark);
    REPORT("Polarisation combination took %s.\n",printed_benchmark);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "polarisation_combination_time",BAD_CAST printed_benchmark);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);

    /* now we search through the observations for the transients that we're looking for */
    STATUS("finding peaks\n");
    START_TIMER;
    peak_finder(&all_transforms,&found_peaks);
    STOP_TIMER;
    TIMER_BENCHMARK;
    PRINT_BENCHMARK(printed_benchmark);
    REPORT("Peak finding took %s.\n",printed_benchmark);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "peak_search_time",BAD_CAST printed_benchmark);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);

    /* and get rid of those peaks that we don't care about */
    STATUS("culling peaks\n");
    START_TIMER;
    cull_peaks(&all_transforms);
    STOP_TIMER;
    TIMER_BENCHMARK;
    PRINT_BENCHMARK(printed_benchmark);
    REPORT("Peak checking took %s.\n",printed_benchmark);
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "peak_verification_time",BAD_CAST printed_benchmark);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);

    /* now separate the data into different lists, one for each source */
    STATUS("separating lists\n");
    n_lists=separate_lists(&all_transforms,&separate_correlations);
    DEBUG("got %d lists\n",n_lists);
    STATUS("making plots\n");
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "PLOT");
    for (j=0;j<n_lists;j++){
      DEBUG("showing list for source %s\n",separate_correlations[j]->sourcename);
      if (!(arguments.plot_type & PLOT_NONE))
	plot_data(&(separate_correlations[j]),arguments.plot_type|PLOT_SEPARATE|PLOT_COMBINED_POLS);
    }
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    /* free our memory */
    FREE(separate_correlations);
    separate_correlations=NULL;
    clear_correlations(&all_transforms);
    FREE(goodfiles[i]);
    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
  }

  /* free all the memory we've allocated for ourselves */
  clear_correlations(&all_correlations);
  clear_correlations(&all_transforms);
  FREE(goodfiles);

  stop_program(0);

  return(0); /* we only have this here so we don't get compiler warnings */
}

void baseline_to_antennae(int baseline_number,int *antenna1,int *antenna2){
  /* this routine takes a baseline number and gets the antenna numbers 
     from it, as baseline_number = (256*antenna1) + antenna2 */
  div_t result;
  result=div(baseline_number,256);
  *antenna1=result.quot;
  *antenna2=result.rem;
}

int want_this_file(int *longest_baseline){
  /* this routine goes through the information in the RPFITS header and
     decides whether the file is the right type for our project.
     it will return 1 if the file is the right type, 0 otherwise */
  int i,retval=1,j,matched;
  char datobs[BUFSIZE],source_name[BUFSIZE],instrument[BUFSIZE],stokesname[BUFSIZE];
  char tmp[BUFSIZE];
  time_t now_time,then_time,cutoff_time;
  double baseline_length,max_baseline_length=0;
  struct tm time_then,time_now;
  
  /* debugging: print the header info */
  DEBUG("number of IFs: %d\n",if_.n_if);
  DEBUG("info         : %s\n",names_.instrument);
  strncpy(instrument,names_.instrument,16);
  for (i=0;i<16;i++){
    if (instrument[i]==' '){
      instrument[i]='\0';
      break;
    }
  }
  instrument[15]='\0';
  strip_trailing_spaces(instrument);
  DEBUG("instrument   : %s\n",instrument);
  strncpy(datobs,names_.datobs,12);
  datobs[10]='\0';
  strip_trailing_spaces(datobs);
  REPORT("This data was obtained at %s, has a cycle time of %d s, and has %d IFs:\n",datobs,
	 param_.intime,if_.n_if);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "date",BAD_CAST datobs);
  snprintf(tmp,BUFSIZE,"%d",param_.intime);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "cycle_time",BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%d",if_.n_if);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "number_ifs",BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%d",su_.n_su);
  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	    BAD_CAST "number_sources",BAD_CAST tmp);
  DEBUG("datobs       : %s\n",datobs);
  DEBUG("feed_type    : %s\n",names_.feed_type);
  for (i=0;i<if_.n_if;i++){
    DEBUG("if_nfreq[%d]  : %d\n",i,if_.if_nfreq[i]);
    REPORT(" IF %d has central frequency %.1f MHz, bandwidth %.1f MHz, %d channels\n",if_.if_num[i],
	   doubles_.if_freq[i]/1E6,doubles_.if_bw[i]/1E6,if_.if_nfreq[i]);
    if (i>0){
      /* have to close the last IF element first */
      XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    }
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "IF");
    snprintf(tmp,BUFSIZE,"%d",if_.if_num[i]);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "number",BAD_CAST tmp);
    snprintf(tmp,BUFSIZE,"%.1f",doubles_.if_freq[i]/1E6);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "central_frequency",BAD_CAST tmp);
    snprintf(tmp,BUFSIZE,"%.1f",doubles_.if_bw[i]/1E6);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "bandwidth",BAD_CAST tmp);
    snprintf(tmp,BUFSIZE,"%d",if_.if_nfreq[i]);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "channels",BAD_CAST tmp);
    DEBUG("if_nstok[%d]  : %d\n",i,if_.if_nstok[i]);
    for (j=0;j<if_.if_nstok[i];j++){
      strncpy(stokesname,names_.if_cstok+i*(if_.if_nstok[i])*2+j*2,2);
      stokesname[2]='\0';
      DEBUG(" if_cstok[%d] : %s\n",j,stokesname);
    }
    DEBUG("if_freq[%d]   : %f\n",i,doubles_.if_freq[i]);
    DEBUG("if_bw[%d]     : %f\n",i,doubles_.if_bw[i]);
    DEBUG("if_num[%d]    : %d\n",i,if_.if_num[i]);
    DEBUG("if_simul[%d]  : %d\n",i,if_.if_simul[i]);
  }
  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
  DEBUG("inttime      : %d\n",param_.intime);
  DEBUG("intbase      : %f\n",param_.intbase);
  DEBUG("# of sources : %d\n",su_.n_su);
  REPORT("This data contains %d sources: ",su_.n_su);
  for (i=0;i<su_.n_su;i++){
    strncpy(source_name,names_.su_name+i*16,16);
    source_name[15]='\0';
    strip_trailing_spaces(source_name);
    DEBUG("source %d name: %s\n",i,source_name);
    if (i>0){
      REPORT(", ");
    }
    REPORT("%s",source_name);
    if (i>0){
      /* have to close the last IF element first */
      XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
    }
    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SOURCE");
    snprintf(tmp,BUFSIZE,"%d",i+1);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "number",BAD_CAST tmp);
    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "name",BAD_CAST source_name);
  }
  REPORT("\n");
  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
  DEBUG("# of ant     : %d\n",anten_.nant);
  for (i=0;i<anten_.nant;i++){
    for (j=0;j<anten_.nant;j++){
      baseline_length=sqrt(pow((doubles_.x[i]-doubles_.x[j]),2)+
			   pow((doubles_.y[i]-doubles_.y[j]),2)+
			   pow((doubles_.z[i]-doubles_.z[j]),2));
      if (baseline_length>max_baseline_length){
	max_baseline_length=baseline_length;
	*longest_baseline=(256*anten_.ant_num[i])+anten_.ant_num[j];
      }
    }
  }
  DEBUG("mx baseline : %f\n",max_baseline_length);
  DEBUG("long base   : %d\n",*longest_baseline);
  
  /* OK, so we want files meeting the following criteria: */
  /* need to be using the ATCA */
  matched=0;
  for (i=0;i<CONSIDER_N_INSTRUMENTS;i++){
    if (strcmp(instrument,consider_instrument[i])==0)
      matched=1;
  }
  if (matched==0){
    DEBUG("incorrect instrument\n");
    retval=0;
  }
  
  /* need to have 2 IFs */
#ifdef MORE_IFS
  if (if_.n_if<MIN_N_IFS)
#else
  if (if_.n_if!=MIN_N_IFS)
#endif
  {
    DEBUG("# of IFS bad\n");
    retval=0;
  }
  
  /* these two IFs needed to be observed simultaneously */
#ifndef NOT_SIMULTANEOUS
  matched=1;
  for (i=1;i<if_.n_if;i++){
    if (if_.if_simul[0]!=if_.if_simul[i])
      matched=0;
  }
  if (matched==0){
    DEBUG("IFs not simultaneous\n");
    retval=0;
  }
#endif
  
  /* the observation date needs to be more than 18 months ago */
  time(&now_time); /* the current time */
  gmtime_r(&now_time,&time_now);
  /* we make a time_t that would be midnight UTC on the day AFTER the observation
     was taken */
  time_then.tm_sec=time_then.tm_min=time_then.tm_hour=0; /* midnight */
  sscanf(datobs,"%4d-%2d-%2d",&(time_then.tm_year),&(time_then.tm_mon),&(time_then.tm_mday));
  time_then.tm_year-=1900; /* years since 1900 */
  time_then.tm_mon-=1; /* months since january */
  time_then.tm_mday+=1; /* the day after the observation */
  then_time=timegm(&time_then);
  /* now we figure out 18 months before now */
  time_now.tm_mon-=18;
  cutoff_time=timegm(&time_now);
  /* now do the comparison */
  if (then_time>cutoff_time){
    DEBUG("observation too recent\n");
    retval=0;
  }
  
  for (i=0;i<if_.n_if;i++){
    /* need to have 32-channel data, ie. continuum observation */
#ifdef LESS_CHANNELS
    if (if_.if_nfreq[i]>(MAX_N_CHANNELS+1)) /* number of frequency channels is one less than that given */
#else
    if (if_.if_nfreq[i]!=(MAX_N_CHANNELS+1))
#endif
    {
      DEBUG("incorrect number of channels\n");
      retval=0;
    }
    /* both IFs need to be at L-band */
    if (doubles_.if_freq[i]>MAX_IF_FREQUENCY){
      DEBUG("- frequencies incorrect\n");
      retval=0;
    }
    /* both IFs need to have a bandwidth greater than 64MHz */
    if (doubles_.if_bw[i]<MIN_BANDWIDTH){
      DEBUG("bandwidth too narrow\n");
      retval=0;
    }
  }

  return (retval);
  
}

int choose_CA_files(const struct dirent *ep){
  /* this routine returns 1 for files which match the pattern *.C*, and 
     returns 0 for all other files */
  int fnmatch_result;
  
  fnmatch_result=fnmatch("*.C*",ep->d_name,0);
  /* don't want C999 files though */
  if (fnmatch_result==0)
    fnmatch_result=!fnmatch("*.C999",ep->d_name,0);
  return (!fnmatch_result);
}

int read_data(correlation_data **all_correlations,int longest_baseline,int *keep_reading){
  /* this routine goes through an RPFITS file and reads in all the 
     data, putting it in a data structure for easy access later */
  int reading_data=-1,current_if=-1,this_jstat=-1,rpfits_result=-1;
  int baseline=-1,scannum=-1,flag=-1,bin=-1,if_no=-1,sourceno=-1;
  int correlations_added=0,add_result;
  float *vis=NULL,*wgt=NULL,ut=0,u=0,v=0,w=0;

  /* now we figure out how much data there is in this scan */
  reading_data=1;
  current_if=-1;
  scannum=0;
  while(reading_data==1){
    this_jstat=JSTAT_READDATA;
    rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,
			    &flag,&bin,&if_no,&sourceno);
    DEBUG("= rpfits_result = %d, baseline = %d, scannum = %d\n",this_jstat,baseline,scannum);
    DEBUG("= flag = %d, if_no = %d, expected if = %d, sourceno = %d\n",
	  flag,if_no,current_if,sourceno);
    if (baseline==-1){
      /* calibration data, don't bother with it */
      /* but this does mean we need to change baselines */
      current_if++;
      if (current_if==if_.n_if){
	current_if=0;
	scannum++;
      }
      /* allocate the vis and wgt arrays to be the correct size */
      if (current_if!=-1){
	if (vis!=NULL){
	  FREE(vis);
	}
	if (wgt!=NULL){
	  FREE(wgt);
	}
	MALLOC(vis,((if_.if_nstok[current_if])*(if_.if_nfreq[current_if]))*2*sizeof(float));
	MALLOC(wgt,((if_.if_nstok[current_if])*(if_.if_nfreq[current_if]))*sizeof(float));
      }
    }
    if ((this_jstat==JSTAT_SUCCESSFUL)||
	(this_jstat==JSTAT_ENDOFSCAN)){
      if ((scannum>=SCAN_BOUNDARY_LIMIT)&&(baseline!=-1)){
	/* add this data to the structure */
	add_result=add_correlation(all_correlations,&vis,&wgt,baseline,ut,u,v,w,flag,bin,if_no,
				   sourceno,if_.if_nfreq[current_if],if_.if_nstok[current_if],
				   longest_baseline,scannum);
	if (add_result==0)
	  correlations_added++;
      }
    } else if (this_jstat==JSTAT_ENDOFFILE){
      reading_data=0;
      *keep_reading=0;
    } else
      reading_data=0;
  }
  FREE(vis);
  FREE(wgt);

  return(correlations_added);
}

void strip_pgplot_device_specifier(char *PGPLOT_device,char *PGPLOT_prefix){
  /* this routine takes a PGPLOT device specifier like file.ps/ps and gets
     rid of the last / and the stuff after it, leaving just the filename
     essentially */
  int i;

  strcpy(PGPLOT_prefix,PGPLOT_device);
  for (i=strlen(PGPLOT_prefix)-1;i>=0;i--){
    if (PGPLOT_prefix[i]=='/'){
      PGPLOT_prefix[i]='\0';
      break;
    }
  }
}

/* how should the peak summary plot be divided up?
   the following settings are to make a 5x6 array of plots */
#define PANEL_MIN_ANT 1
#define PANEL_MAX_ANT 6
#define PANEL_X_COLS 5
#define PANEL_Y_ROWS 6
#define PANEL_MIN_IF 1
#define PANEL_MAX_IF 2

#define PGPLOT_MEDIUM_SCREEN -1
#define PGPLOT_MEDIUM_PAPER -2

int plot_data(correlation_data **all_correlations,int plot_type){
  /* this routine plots the data it gets sent in one of the ways that it can */
  int k,l,m,n_times,current_waterfall_channel,i,j;
  int current_waterfall_index,npoints,*symbols=NULL,ant1,ant2;
  int xpanel,ypanel,made_plot=0,tmpchannel;
  correlation_data *cycle_data=NULL,*second_cycle=NULL,*cycle_hold=NULL;
  correlation_data *another_cycle=NULL;
  float minreal,maxreal,minimag,maximag,*channels=NULL,***spectrum=NULL;
  float *waterfall_array=NULL,tr[6];
  float plot_minvp_x,plot_minvp_y,plot_maxvp_x,plot_maxvp_y,plot_nearby_height;
  float plot_main_buffer,plot_nearby_buffer,plot_nearby_width;
  float straightline_x[2],straightline_y[2],*obs_times=NULL,*avg_amps=NULL;
  float *top_errs=NULL,*bottom_errs=NULL,plot_row_height,plot_col_width;
  float baseline_length,x_string_length,y_string_length,old_height;
  char title[BUFSIZE],PGPLOT_device[BUFSIZE],timestring[BUFSIZE],tmpfile[BUFSIZE];
  char baseline_string[BUFSIZE],tmp[BUFSIZE];
  time_t time_difference;
  double mintime,maxtime;

  if (*all_correlations==NULL)
    /* nothing to plot */
    return(0);

  /* interpret the PGPLOT device name */
  /* look for the symbols: %S = source name 
                           %F = RPFITS file name */
  PGPLOT_device[0]='\0';
  for (i=0;i<strlen(arguments.PGPLOT_device);i++){
    if (arguments.PGPLOT_device[i]!='%'){
      j=strlen(PGPLOT_device);
      PGPLOT_device[j]=arguments.PGPLOT_device[i];
      PGPLOT_device[j+1]='\0';
    } else {
      i++;
      if (arguments.PGPLOT_device[i]=='S'){
	/* need to replace this with the source name */
	strcat(PGPLOT_device,(*all_correlations)->sourcename);
      } else if (arguments.PGPLOT_device[i]=='F'){
	/* need to replace this with the RPFITS file name */
	strip_path(names_.file,tmpfile);
	strcat(PGPLOT_device,tmpfile);
      }
    }
  }
  
  DEBUG("using PGPLOT device %s\n",PGPLOT_device);
  if (strcmp(opened_pgplot_device_name,PGPLOT_device)!=0){
    if (strlen(opened_pgplot_device_name)!=0){
      cpgclos();
      opened_pgplot_device_name[0]='\0';
    }
    opened_pgplot_device_id=cpgopen(PGPLOT_device);
    strcpy(opened_pgplot_device_name,PGPLOT_device);
    cpgask(arguments.PGPLOT_ask);
  } else {
    cpgslct(opened_pgplot_device_id);
  }

  if (plot_type & PLOT_INDIVIDUAL_SPECTRA){
    /* now go through the data correlation by correlation */
    cycle_data=*all_correlations;
    while(cycle_data!=NULL){
      if (cycle_data->islongest_baseline==1){
	minreal=minimag=1E6;
	maxreal=maximag=-1E6;
	MALLOC(channels,cycle_data->if_nfreq*sizeof(float));
	if (plot_type & PLOT_COMBINED_POLS){
	  MALLOC(spectrum,1*sizeof(float**));
	} else {
	  MALLOC(spectrum,cycle_data->if_nstok*sizeof(float**));
	}
	for (k=0;k<cycle_data->if_nstok;k++){
	  if ((plot_type & PLOT_COMBINED_POLS)&&(k==0)){
	    MALLOC(spectrum[k],1*sizeof(float*));
	  } else if (!(plot_type & PLOT_COMBINED_POLS)) {
	    MALLOC(spectrum[k],2*sizeof(float*));
	  }
	  for (l=0;l<2;l++){
	    if ((!(plot_type & PLOT_COMBINED_POLS))||((plot_type & PLOT_COMBINED_POLS)&&(k==0)&&(l==0))){
	      MALLOC(spectrum[k][l],cycle_data->if_nfreq*sizeof(float));
	    }
	    for (m=0;m<cycle_data->if_nfreq;m++){
	      if ((plot_type & PLOT_COMBINED_POLS)&&(k==0)&&(l==0))
		spectrum[k][l][m]=cycle_data->pols_combined[m];
	      else if (!(plot_type & PLOT_COMBINED_POLS))
		spectrum[k][l][m]=cycle_data->vis[m*2*cycle_data->if_nstok+k*2+l];
	      if ((!(plot_type & PLOT_COMBINED_POLS))||((plot_type & PLOT_COMBINED_POLS)&&(k==0)&&(l==0))){
		if ((l==0)&&(spectrum[k][l][m]<minreal))
		  minreal=spectrum[k][l][m];
		if ((l==0)&&(spectrum[k][l][m]>maxreal))
		  maxreal=spectrum[k][l][m];
		if ((l==1)&&(spectrum[k][l][m]<minimag))
		  minimag=spectrum[k][l][m];
		if ((l==1)&&(spectrum[k][l][m]>maximag))
		  maximag=spectrum[k][l][m];
		channels[m]=m;
	      }
	    }
	  }
	}
	if (plot_type & PLOT_COMBINED_POLS)
	  snprintf(title,BUFSIZE,"%s baseline %d scannum %d flag %d if_no %d Stokes %s",
		   cycle_data->sourcename,cycle_data->baseline,cycle_data->scannum,
		   cycle_data->flag,cycle_data->if_no,cycle_data->combined_pol_name);
	else
	  snprintf(title,BUFSIZE,"%s baseline %d scannum %d flag %d if_no %d",
		   cycle_data->sourcename,cycle_data->baseline,cycle_data->scannum,
		   cycle_data->flag,cycle_data->if_no);
	cpgpage();
	made_plot=1;
	if (plot_type & PLOT_COMBINED_POLS)
	  cpgsvp(0.12,0.98,0.12,0.94);
	else
	  cpgsvp(0.12,0.98,0.55,0.88);
	cpgswin(0,cycle_data->if_nfreq,minreal,maxreal);
	if (plot_type & PLOT_COMBINED_POLS)
	  cpgbox("BCNTS",0,0,"BCNTSV",0,0);
	else
	  cpgbox("BCTS",0,0,"BCNTSV",0,0);
	cpgline(cycle_data->if_nfreq,channels,spectrum[0][0]);
	if (plot_type & PLOT_COMBINED_POLS)
	  cpglab("channel","amplitude",title);
	else
	  cpglab("real channel","amplitude",title);
	if (!(plot_type & PLOT_COMBINED_POLS)){
	  cpgsvp(0.12,0.98,0.12,0.45);
	  cpgswin(0,cycle_data->if_nfreq,minimag,maximag);
	  cpgbox("BCTS",0,0,"BCNTSV",0,0);
	  cpgline(cycle_data->if_nfreq,channels,spectrum[0][1]);
	  cpglab("imaginary channel","amplitude","");
	}
	/* free some memory */
	if (spectrum!=NULL){
	  for (k=0;k<cycle_data->if_nstok;k++){
	    for (l=0;l<2;l++){
	      if ((!(plot_type & PLOT_COMBINED_POLS))||((plot_type & PLOT_COMBINED_POLS)&&(k==0)&&(l==0)))
		FREE(spectrum[k][l]);
	    }
	    if ((!(plot_type & PLOT_COMBINED_POLS))||((plot_type & PLOT_COMBINED_POLS)&&(k==0)))
	      FREE(spectrum[k]);
	  }
	  FREE(spectrum);
	}
	if (channels!=NULL){
	  FREE(channels);
	}
      }
      if (plot_type & PLOT_ALL)
	cycle_data=cycle_data->next;
      else if (plot_type & PLOT_SEPARATE)
	cycle_data=cycle_data->separate_next;
    }
  } else if (plot_type & PLOT_WATERFALL){
    /* we need to have an array which is number_of_channels x number_of_timesteps 
       in size, so we need to determine the number of timesteps */
    minreal=1E6;
    maxreal=-1E6;
    if (plot_type & PLOT_ALL)
      for (cycle_data=*all_correlations,n_times=0;cycle_data!=NULL;
	   cycle_data=cycle_data->next,n_times++);
    else if (plot_type & PLOT_SEPARATE)
      for (cycle_data=*all_correlations,n_times=0;cycle_data!=NULL;
	   cycle_data=cycle_data->separate_next,n_times++);

    MALLOC(waterfall_array,(n_times* (*all_correlations)->if_nfreq)*sizeof(float));
    /* now go through the data and assign values to the array */
    current_waterfall_index=0;
    cycle_data=*all_correlations;
    while(cycle_data!=NULL){
      for (current_waterfall_channel=0;current_waterfall_channel<cycle_data->if_nfreq;
	   current_waterfall_channel++){
	if (plot_type & PLOT_COMBINED_POLS)
	  waterfall_array[current_waterfall_index]=cycle_data->pols_combined[current_waterfall_channel];
	else
	  waterfall_array[current_waterfall_index]=cycle_data->vis[current_waterfall_channel*2*cycle_data->if_nstok];
 	if (waterfall_array[current_waterfall_index]<minreal)
 	  minreal=waterfall_array[current_waterfall_index];
 	if (waterfall_array[current_waterfall_index]>maxreal)
 	  maxreal=waterfall_array[current_waterfall_index];
 	current_waterfall_index++;
      }
      if (plot_type & PLOT_ALL)
	cycle_data=cycle_data->next;
      else if (plot_type & PLOT_SEPARATE) 
	cycle_data=cycle_data->separate_next; 
    }
    tr[0]=tr[2]=tr[3]=tr[4]=0.0;
    tr[1]=tr[5]=1.0;

    /* set up a colourmap first */
    cpgswin(0,(*all_correlations)->if_nfreq+1,0,n_times+1);
    cpgscir(0,2000);
    cpgsitf(0);
    palett(PGPLOT_CT_GREYSCALE,1.0,0.5);
    cpgimag(waterfall_array,(*all_correlations)->if_nfreq,n_times,1,(*all_correlations)->if_nfreq,
 	    1,n_times,minreal,maxreal,tr);
    made_plot=1;
    FREE(waterfall_array);
  } else if (plot_type & PLOT_PEAKS){
    /* we plot spectra that have large peaks in the lag spectra */
    /* go through the data and find a spectrum that displays a peak */
    cycle_data=*all_correlations;
    while(cycle_data!=NULL){
      if (cycle_data->spectrum_info.npeaks_found>0){
	/* we assign the boundaries of the viewport */
	plot_minvp_x=0.05;
	plot_minvp_y=0.02;
	plot_maxvp_x=0.93;
	plot_maxvp_y=0.93;
	plot_nearby_height=0.17;
	plot_nearby_buffer=0.08;
	plot_main_buffer=0.15;
	cpgpage();
	made_plot=1;
	if ((cycle_data->islongest_baseline==1)&&
	    (cycle_data->if_no==1)){
	  /* put the details of the plot in the report */
	  if (output_medium(PGPLOT_device)==PGPLOT_MEDIUM_PAPER){
	    strip_pgplot_device_specifier(PGPLOT_device,tmpfile);
	  } else {
	    strcpy(tmpfile,"screen");
	  }
	  REPORT("Plotting peak P%04d to %s\n",
		 cycle_data->spectrum_info.peak_number[0],tmpfile);
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "PEAK");
	  snprintf(tmp,BUFSIZE,"%04d",cycle_data->spectrum_info.peak_number[0]);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "number",BAD_CAST tmp);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "filename",BAD_CAST tmpfile);
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	}

	/* make the channels array */
	MALLOC(channels,cycle_data->if_nfreq*sizeof(float));
	for (i=0;i<cycle_data->if_nfreq;i++)
	  channels[i]=i;
	/* first we plot the spectra of the data before and after it */
	/* calculate the width of each plot */
	plot_nearby_width=((plot_maxvp_x-plot_minvp_x
			    -((arguments.number_of_nearby_spectra-1)*plot_nearby_buffer))/
			   arguments.number_of_nearby_spectra);
 	DEBUG("plot_nearby_width = %f\n",plot_nearby_width);
	/* we'll use the same axis ranges for all plots, and we'll take that from
	   this spectrum */
	minreal=1E6;
	maxreal=-1E6;
	for (i=0;i<cycle_data->if_nfreq;i++){
	  if (cycle_data->pols_combined[i]<minreal)
	    minreal=cycle_data->pols_combined[i];
	  if (cycle_data->pols_combined[i]>maxreal)
	    maxreal=cycle_data->pols_combined[i];
	}
	straightline_x[0]=0;
	straightline_x[1]=cycle_data->if_nfreq;
	/* before it first, keeping to the same baseline */
	for (second_cycle=cycle_data->samebaseline_prev,i=0;
	     ((second_cycle!=NULL)&&(i<arguments.number_of_nearby_spectra));
	     second_cycle=second_cycle->samebaseline_prev,i++){
	  /* set where the plot will show up */
 	  DEBUG("vp coords = %f, %f -> %f, %f\n",
		plot_maxvp_x-(i+1)*plot_nearby_width-i*plot_nearby_buffer,
		plot_maxvp_y-plot_nearby_height,
		plot_maxvp_x-i*plot_nearby_width-i*plot_nearby_buffer,
		plot_maxvp_y);
	  cpgsvp(plot_maxvp_x-(i+1)*plot_nearby_width-i*plot_nearby_buffer,
		 plot_maxvp_x-i*plot_nearby_width-i*plot_nearby_buffer,
		 plot_maxvp_y-plot_nearby_height,plot_maxvp_y);
	  /* set the axis range */
	  cpgswin(0,second_cycle->if_nfreq,minreal,maxreal);
	  /* draw the axes */
	  cpgsci(PGSCI_INVERSE);
	  cpgsls(PGSLS_FULL);
	  cpgbox("BCTS",0,0,"BCTS",0,0);
	  /* draw the spectrum */
	  cpgline(second_cycle->if_nfreq,channels,second_cycle->pols_combined);
	  /* plot the average power */
	  straightline_y[0]=straightline_y[1]=second_cycle->spectrum_info.average_amplitude;
	  cpgsci(PGSCI_RED);
	  cpgsls(PGSLS_FULL);
	  cpgline(2,straightline_x,straightline_y);
	  /* and a standard deviation above it */
	  straightline_y[0]=straightline_y[1]=second_cycle->spectrum_info.average_amplitude
	    +second_cycle->spectrum_info.rms_noise;
	  cpgsci(PGSCI_GREEN);
	  cpgsls(PGSLS_DASHED);
	  cpgline(2,straightline_x,straightline_y);
	  /* now set the title to be how long this spectra was taken before the main one */
	  time_difference=second_cycle->observation_time-cycle_data->observation_time;
 	  DEBUG("sc obs time = %ld, cd obs time = %ld diff = %ld\n",second_cycle->observation_time,
 		 cycle_data->observation_time,time_difference);
	  snprintf(title,BUFSIZE,"%+ld seconds",time_difference);
	  cpgsci(PGSCI_INVERSE);
	  cpgsls(PGSLS_FULL);
	  cpglab("","",title);
	}
	/* now after it */
	for (second_cycle=cycle_data->samebaseline_next,i=0;
	     ((second_cycle!=NULL)&&(i<arguments.number_of_nearby_spectra));
	     second_cycle=second_cycle->samebaseline_next,i++){
	  /* set where the plot will show up */
	  cpgsvp(plot_maxvp_x-(i+1)*plot_nearby_width-i*plot_nearby_buffer,
		 plot_maxvp_x-i*plot_nearby_width-i*plot_nearby_buffer,
		 plot_minvp_y,plot_minvp_y+plot_nearby_height);
	  /* set the axis range */
	  cpgswin(0,second_cycle->if_nfreq,minreal,maxreal);
	  /* draw the axes */
	  cpgsci(PGSCI_INVERSE);
	  cpgsls(PGSLS_FULL);
	  cpgbox("BCTS",0,0,"BCTS",0,0);
	  /* draw the spectrum */
	  cpgline(second_cycle->if_nfreq,channels,second_cycle->pols_combined);
	  /* plot the average power */
	  straightline_y[0]=straightline_y[1]=second_cycle->spectrum_info.average_amplitude;
	  cpgsci(PGSCI_RED);
	  cpgsls(PGSLS_FULL);
	  cpgline(2,straightline_x,straightline_y);
	  /* and a standard deviation above it */
	  straightline_y[0]=straightline_y[1]=second_cycle->spectrum_info.average_amplitude
	    +second_cycle->spectrum_info.rms_noise;
	  cpgsci(PGSCI_GREEN);
	  cpgsls(PGSLS_DASHED);
	  cpgline(2,straightline_x,straightline_y);
	  /* now set the title to be how long this spectra was taken before the main one */
	  time_difference=second_cycle->observation_time-cycle_data->observation_time;
	  snprintf(title,BUFSIZE,"%+ld seconds",time_difference);
	  cpgsci(PGSCI_INVERSE);
	  cpgsls(PGSLS_FULL);
	  cpglab("","",title);
	}
	/* now the main plot, with the peak identified */
	/* set where the plot will show up */
	cpgsvp(plot_maxvp_x-0.5*(plot_maxvp_x-plot_minvp_x),plot_maxvp_x,
	       plot_minvp_y+plot_nearby_height+plot_main_buffer,
	       plot_maxvp_y-plot_nearby_height-plot_main_buffer);
	/* set the axis range */
	cpgswin(0,cycle_data->if_nfreq,minreal,maxreal);
	/* draw the axes */
	cpgsci(PGSCI_INVERSE);
	cpgsls(PGSLS_FULL);
	cpgbox("BCNTS",0,0,"BCNTS",0,0);
	/* draw the spectrum */
	cpgline(cycle_data->if_nfreq,channels,cycle_data->pols_combined);
	/* plot the average power that the spectrum says it has */
	straightline_y[0]=straightline_y[1]=cycle_data->spectrum_info.average_amplitude;
	cpgsci(PGSCI_RED);
	cpgsls(PGSLS_FULL);
	cpgline(2,straightline_x,straightline_y);
	/* also plot the average power that was estimated from the nearby spectra */
	straightline_y[0]=straightline_y[1]=cycle_data->spectrum_info.estimated_average_amplitude;
	cpgsci(PGSCI_RED);
	cpgsls(PGSLS_DOTDASH);
	cpgline(2,straightline_x,straightline_y);
	/* and a standard deviation above the latter line */
	straightline_y[0]=straightline_y[1]=cycle_data->spectrum_info.estimated_average_amplitude+
	  cycle_data->spectrum_info.estimated_rms_noise;
	cpgsci(PGSCI_GREEN);
	cpgsls(PGSLS_DASHED);
	cpgline(2,straightline_x,straightline_y);
	/* now draw vertical lines along each channel where there is an identified peak */
	straightline_y[0]=minreal;
	straightline_y[1]=maxreal;
	if (output_medium(PGPLOT_device)==PGPLOT_MEDIUM_PAPER)
	  /* 	if (strncmp(PGPLOT_device+strlen(PGPLOT_device)-2,"ps",2)==0) */
	  /* this is a paper output, use a dark colour */
	  cpgsci(PGSCI_BLUE);
	else
	  /* this is probably a screen output, use a light colour */
	  cpgsci(PGSCI_YELLOW);
	cpgsls(PGSLS_DASHED);
	for (i=0;i<cycle_data->spectrum_info.npeaks_found;i++){
	  straightline_x[0]=straightline_x[1]=cycle_data->spectrum_info.peak_channels[i];
	  if (cycle_data->spectrum_info.peak_present_both_pols[i]==0)
	    cpgsci(PGSCI_RED);
	  cpgline(2,straightline_x,straightline_y);
	}
	/* now set the title to be the usual information */
	baseline_to_antennae(cycle_data->baseline,&ant1,&ant2);
	strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(cycle_data->observation_time)));
	strip_path(names_.file,tmpfile);
	snprintf(title,BUFSIZE,"%s P%04d %s ants (%d, %d) @ %s if_no %d Stokes %s",
		 tmpfile,cycle_data->spectrum_info.peak_number[0],cycle_data->sourcename,
		 ant1,ant2,timestring,cycle_data->if_no,cycle_data->combined_pol_name);
	/* and generate the axis labels now */
	cpgsci(PGSCI_INVERSE);
	cpglab("channel","amplitude","");
	/* centre the title a little better */
	cpgsvp(plot_minvp_x,plot_maxvp_x,plot_minvp_y+plot_nearby_height+plot_main_buffer,
	       plot_maxvp_y-plot_nearby_height-plot_main_buffer);
	cpgswin(0,100,0,100);
	cpgptxt(50,130,0.0,0.5,title);
	/* now we also want to make a plot showing how the average amplitudes of the
	   nearby spectra were fitted to get the estimate of the average amplitude of
	   the main spectrum */
	/* set where the plot will show up */
	cpgsvp(plot_minvp_x,plot_maxvp_x-0.5*(plot_maxvp_x-plot_minvp_x)-plot_main_buffer,
	       plot_minvp_y+plot_nearby_height+plot_main_buffer,
	       plot_maxvp_y-plot_nearby_height-plot_main_buffer);
	/* determine the axis range and make the arrays to plot */
	minreal=maxreal=cycle_data->spectrum_info.estimated_average_amplitude;
	mintime=maxtime=cycle_data->observation_time-cycle_data->observation_time;
 	DEBUG("min, max time = %f, %f min, max real = %f, %f\n",
	      mintime,maxtime,minreal,maxreal);
	npoints=1;
	REALLOC(obs_times,npoints*sizeof(float));
	REALLOC(avg_amps,npoints*sizeof(float));
	REALLOC(top_errs,npoints*sizeof(float));
	REALLOC(bottom_errs,npoints*sizeof(float));
	REALLOC(symbols,npoints*sizeof(int));
	obs_times[npoints-1]=cycle_data->observation_time-cycle_data->observation_time;
	avg_amps[npoints-1]=cycle_data->spectrum_info.estimated_average_amplitude;
	top_errs[npoints-1]=cycle_data->spectrum_info.estimated_average_amplitude+
	  cycle_data->spectrum_info.estimated_average_amplitude_error;
	bottom_errs[npoints-1]=cycle_data->spectrum_info.estimated_average_amplitude-
	  cycle_data->spectrum_info.estimated_average_amplitude_error;
	symbols[npoints-1]=-1;
	for (second_cycle=cycle_data->samebaseline_prev,i=0;
	     ((second_cycle!=NULL)&&(i<arguments.number_of_nearby_spectra));
	     second_cycle=second_cycle->samebaseline_prev,i++){
	  npoints++;
	  REALLOC(obs_times,npoints*sizeof(float));
	  REALLOC(avg_amps,npoints*sizeof(float));
	  REALLOC(top_errs,npoints*sizeof(float));
	  REALLOC(bottom_errs,npoints*sizeof(float));
	  REALLOC(symbols,npoints*sizeof(int));
	  obs_times[npoints-1]=second_cycle->observation_time-cycle_data->observation_time;
 	  DEBUG("earlier obs time = %ld\n",second_cycle->observation_time-cycle_data->observation_time);
	  avg_amps[npoints-1]=second_cycle->spectrum_info.average_amplitude;
	  top_errs[npoints-1]=second_cycle->spectrum_info.average_amplitude+
	    second_cycle->spectrum_info.rms_noise;
	  bottom_errs[npoints-1]=second_cycle->spectrum_info.average_amplitude-
	    second_cycle->spectrum_info.rms_noise;
	  symbols[npoints-1]=0;
	  if ((second_cycle->spectrum_info.average_amplitude-
	       second_cycle->spectrum_info.rms_noise)<minreal)
	    minreal=second_cycle->spectrum_info.average_amplitude-
	      second_cycle->spectrum_info.rms_noise;
	  if ((second_cycle->spectrum_info.average_amplitude+
	       second_cycle->spectrum_info.rms_noise)>maxreal)
	    maxreal=second_cycle->spectrum_info.average_amplitude+
	      second_cycle->spectrum_info.rms_noise;
	  if ((second_cycle->observation_time-cycle_data->observation_time)<mintime)
	    mintime=second_cycle->observation_time-cycle_data->observation_time;
	  if ((second_cycle->observation_time-cycle_data->observation_time)>maxtime)
	    maxtime=second_cycle->observation_time-cycle_data->observation_time;
	}
	for (second_cycle=cycle_data->samebaseline_next,i=0;
	     ((second_cycle!=NULL)&&(i<arguments.number_of_nearby_spectra));
	     second_cycle=second_cycle->samebaseline_next,i++){
	  npoints++;
	  REALLOC(obs_times,npoints*sizeof(float));
	  REALLOC(avg_amps,npoints*sizeof(float));
	  REALLOC(top_errs,npoints*sizeof(float));
	  REALLOC(bottom_errs,npoints*sizeof(float));
	  REALLOC(symbols,npoints*sizeof(int));
	  obs_times[npoints-1]=second_cycle->observation_time-cycle_data->observation_time;
 	  DEBUG("later obs time = %ld\n",second_cycle->observation_time-cycle_data->observation_time);
	  avg_amps[npoints-1]=second_cycle->spectrum_info.average_amplitude;
	  top_errs[npoints-1]=second_cycle->spectrum_info.average_amplitude+
	    second_cycle->spectrum_info.rms_noise;
	  bottom_errs[npoints-1]=second_cycle->spectrum_info.average_amplitude-
	    second_cycle->spectrum_info.rms_noise;
	  symbols[npoints-1]=0;
	  if ((second_cycle->spectrum_info.average_amplitude-
	       second_cycle->spectrum_info.rms_noise)<minreal)
	    minreal=second_cycle->spectrum_info.average_amplitude-
	      second_cycle->spectrum_info.rms_noise;
	  if ((second_cycle->spectrum_info.average_amplitude+
	       second_cycle->spectrum_info.rms_noise)>maxreal)
	    maxreal=second_cycle->spectrum_info.average_amplitude+
	      second_cycle->spectrum_info.rms_noise;
	  if ((second_cycle->observation_time-cycle_data->observation_time)<mintime)
	    mintime=second_cycle->observation_time-cycle_data->observation_time;
	  if ((second_cycle->observation_time-cycle_data->observation_time)>maxtime)
	    maxtime=second_cycle->observation_time-cycle_data->observation_time;
	}
 	DEBUG("min, max time = %f, %f min, max real = %f, %f\n",
	      mintime,maxtime,minreal,maxreal);
	cpgswin(mintime,maxtime,minreal,maxreal);
	/* draw the axes */
	cpgsci(PGSCI_INVERSE);
	cpgsls(PGSLS_FULL);
	cpgbox("BCNTS",0,0,"BCNTS",0,0);
	/* draw the points */
	cpgpnts(npoints,obs_times,avg_amps,symbols,npoints);
	/* and the error bars */
	cpgerry(npoints,obs_times,top_errs,bottom_errs,1.0);
	/* now draw the line of best fit */
	straightline_x[0]=mintime;
	straightline_x[1]=maxtime;
	straightline_y[0]=cycle_data->spectrum_info.fit_c0+
	  cycle_data->spectrum_info.fit_c1*straightline_x[0];
	straightline_y[1]=cycle_data->spectrum_info.fit_c0+
	  cycle_data->spectrum_info.fit_c1*straightline_x[1];
 	DEBUG("x[0]=%f x[1]=%f y[0]=%f y[1]=%f\n",straightline_x[0],straightline_x[1],
	      straightline_y[0],straightline_y[1]);
	cpgsci(PGSCI_RED);
	cpgsls(PGSLS_FULL);
	cpgline(2,straightline_x,straightline_y);
	/* now the axis labels */
	cpgsci(PGSCI_INVERSE);
	cpgsci(PGSLS_FULL);
	cpglab("time","average amplitude","");
	/* do we need to make a peak summary page? */
	if ((cycle_data->islongest_baseline)&&
	    (cycle_data->if_no==1)){
	  /* a peak summary page comes after the longest baseline IF 1
	     spectrum page; this page shows all baselines and IFs that are
	     unflagged for that observation time, and for the time immediately
	     before and after
	     the baselines are plotted like this: 
	     12 13 14 15 16
	     23 24 25 26 34  IF1
	     35 36 45 46 56
	     12 13 14 15 16
	     23 24 25 26 34  IF2
	     35 36 45 46 56
	  */
	  /* we assign the boundaries of the viewport */
	  plot_minvp_x=0.05;
	  plot_minvp_y=0.05;
	  plot_maxvp_x=0.98;
	  plot_maxvp_y=0.93;
	  /* set up a loop for the time range */
	  for (j=0;j<3;j++){
	    cpgpage();
	    /* the first thing to do is determine the axis limits, as every panel 
	       will have the same limits */
	    minreal=maxreal=cycle_data->pols_combined[0];
	    for (second_cycle=*all_correlations;second_cycle!=NULL;
		 second_cycle=second_cycle->next){
	      /* check that its the same observation */
	      if ((second_cycle->observation_time==cycle_data->observation_time)&&
		  (strcmp(second_cycle->sourcename,cycle_data->sourcename)==0)){
		/* do we need to move backwards/forwards in time */
		cycle_hold=second_cycle;
		if (j==1){
		  /* we'll move backwards */
		  cycle_hold=second_cycle;
		  second_cycle=second_cycle->samebaseline_prev;
		} else if (j==2) {
		  /* we'll move forwards */
		  cycle_hold=second_cycle;
		  second_cycle=second_cycle->samebaseline_next;
		}
		for (i=0;i<second_cycle->if_nfreq;i++){
		  if (second_cycle->pols_combined[i]<minreal)
		    minreal=second_cycle->pols_combined[i];
		  if (second_cycle->pols_combined[i]>maxreal)
		    maxreal=second_cycle->pols_combined[i];
		}
		second_cycle=cycle_hold;
	      }
	    }
	    /* some plot parameters */
	    plot_row_height=(plot_maxvp_y-plot_minvp_y)/PANEL_Y_ROWS;
	    plot_col_width=(plot_maxvp_x-plot_minvp_x)/PANEL_X_COLS;
	    /* now we go through this same data and make the plots */
	    for (second_cycle=*all_correlations;second_cycle!=NULL;
		 second_cycle=second_cycle->next){
	      /* check that its the same observation */
	      if ((second_cycle->observation_time==cycle_data->observation_time)&&
		  (strcmp(second_cycle->sourcename,cycle_data->sourcename)==0)){
		/* do we need to move backwards/forwards in time */
		cycle_hold=second_cycle;
		if (j==1){
		  /* we'll move backwards */
		  cycle_hold=second_cycle;
		  second_cycle=second_cycle->samebaseline_prev;
		} else if (j==2) {
		  /* we'll move forwards */
		  cycle_hold=second_cycle;
		  second_cycle=second_cycle->samebaseline_next;
		}
		/* which panel are we going to plot this spectrum in? */
		which_panel(second_cycle->baseline,second_cycle->if_no,&xpanel,&ypanel);
		DEBUG("baseline is %d, xpanel = %d ypanel = %d\n",
		      second_cycle->baseline,xpanel,ypanel);
		/* now set the viewport */
		cpgsvp(plot_minvp_x+xpanel*plot_col_width,
		       plot_minvp_x+(xpanel+1)*plot_col_width,
		       plot_minvp_y+ypanel*plot_row_height,
		       plot_minvp_y+(ypanel+1)*plot_row_height);
		/* and set the axis ranges */
		cpgswin(0,second_cycle->if_nfreq,minreal,maxreal);
		/* draw the axes */
		cpgsci(PGSCI_INVERSE);
		cpgsls(PGSLS_FULL);
		if ((xpanel==0)&&(ypanel==0)){
		  /* we want to label the amplitude and channel axes */
		  cpgbox("BCNTS",0,0,"BCNTS",0,0);
		} else {
		  cpgbox("BCTS",0,0,"BCTS",0,0);
		}
		/* draw the spectrum */
		cpgline(second_cycle->if_nfreq,channels,second_cycle->pols_combined);
		/* draw the average power */
		straightline_x[0]=0;
		straightline_x[1]=second_cycle->if_nfreq;
		straightline_y[0]=straightline_y[1]=
		  second_cycle->spectrum_info.estimated_average_amplitude;
		cpgsci(PGSCI_RED);
		cpgsls(PGSLS_DOTDASH);
		cpgline(2,straightline_x,straightline_y);
		/* and a standard deviation above the line */
		straightline_y[0]=straightline_y[1]=
		  second_cycle->spectrum_info.estimated_average_amplitude+
		  second_cycle->spectrum_info.estimated_rms_noise;
		cpgsci(PGSCI_GREEN);
		cpgsls(PGSLS_DASHED);
		cpgline(2,straightline_x,straightline_y);
		/* and a vertical line along each peak channel */
		straightline_y[0]=minreal;
		straightline_y[1]=maxreal;
		if (output_medium(PGPLOT_device)==PGPLOT_MEDIUM_PAPER)
/* 		if (strncmp(PGPLOT_device+strlen(PGPLOT_device)-2,"ps",2)==0) */
		  /* this is a paper output, use a dark colour */
		  cpgsci(PGSCI_BLUE);
		else
		  /* this is probably a screen output, use a light colour */
		  cpgsci(PGSCI_YELLOW);
		cpgsls(PGSLS_DASHED);
		for (i=0;i<second_cycle->spectrum_info.npeaks_found;i++){
		  straightline_x[0]=straightline_x[1]=second_cycle->spectrum_info.peak_channels[i];
		  /* change the colour to red if the peak wasn't detected in
		     both polarisations */
		  if (second_cycle->spectrum_info.peak_present_both_pols[i]==0)
		    cpgsci(PGSCI_RED);
		  cpgline(2,straightline_x,straightline_y);
		}
		/* we should plot where the peaks should have been in 
		   another colour */
		/* find the longest baseline */
		another_cycle=second_cycle;
		while((another_cycle->observation_time==second_cycle->observation_time)&&
		      (another_cycle->islongest_baseline==0)){
		  another_cycle=another_cycle->prev;
		}
		if ((another_cycle->observation_time!=second_cycle->observation_time)||
		    (another_cycle->islongest_baseline==0)){
		  /* we didn't find the longest baseline going backwards so instead
		     we go forwards */
		  another_cycle=second_cycle;
		  while((another_cycle->observation_time==second_cycle->observation_time)&&
			(another_cycle->islongest_baseline==0)){
		    another_cycle=another_cycle->next;
		  }
		}
		/* we now assume we have found the correct long baseline! 
		   and we also only look at the first peak on this baseline */
		if (another_cycle->spectrum_info.peak_channels!=NULL){
		  lag_scaler(another_cycle->if_bandwidth,another_cycle->spectrum_info.peak_channels[0],
			     another_cycle->u,another_cycle->v,another_cycle->w,
			     second_cycle->u,second_cycle->v,second_cycle->w,
			     NULL,NULL,&tmpchannel);
		  straightline_x[1]=straightline_x[0]=tmpchannel;
		  straightline_y[0]=minreal;
		  straightline_y[1]=maxreal;
		  if (output_medium(PGPLOT_device)==PGPLOT_MEDIUM_PAPER)
		    cpgsci(PGSCI_ORANGE);
		  else
		    cpgsci(PGSCI_ORANGE);
		  cpgsls(PGSLS_DOTTED);
		  cpgslw(3);
		  cpgline(2,straightline_x,straightline_y);
		  cpgslw(1);
		}
		/* now put the baseline length in the panel */
		baseline_length=sqrtf(powf(second_cycle->u,2)+powf(second_cycle->v,2)+
				      powf(second_cycle->w,2));
		baseline_to_antennae(second_cycle->baseline,&ant1,&ant2);
		snprintf(baseline_string,BUFSIZE,"(%d,%d) %.2f m",ant1,ant2,baseline_length);
		/* get the length of the baseline string so we know how to position it */
		cpgqch(&old_height);
		cpgsch(0.6);
		cpglen(PGLEN_WORLD,baseline_string,&x_string_length,&y_string_length);
		/* now we place the string at the top right of the viewport */
		/* choose the colour based on whether the peak was found in both of the
		   individual polarisations */
		cpgsci(PGSCI_INVERSE);
		cpgtext((float)second_cycle->if_nfreq-1.2*x_string_length,
			maxreal-(maxreal-minreal)/8.0,baseline_string);
		cpgsch(old_height);
		second_cycle=cycle_hold;
	      }
	    }
	    /* now we put a title on the plot */
	    if (j==0){
	      strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(cycle_data->observation_time)));
	    } else if (j==1){
	      strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(cycle_data->samebaseline_prev->observation_time)));
	    } else if (j==2){
	      strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(cycle_data->samebaseline_next->observation_time)));
	    }
	    strip_path(names_.file,tmpfile);
	    snprintf(title,BUFSIZE,"%s P%04d %s @ %s Stokes %s",
		     tmpfile,cycle_data->spectrum_info.peak_number[0],cycle_data->sourcename,
		     timestring,cycle_data->combined_pol_name);
	    cpgsvp(plot_minvp_x,plot_maxvp_x,plot_minvp_y,plot_maxvp_y);
	    cpgswin(0,100,0,100);
	    cpgsci(PGSCI_INVERSE);
	    cpglab("","",title);
	  }
	}
	/* free the memory we've allocated */
	FREE(obs_times);
	FREE(avg_amps);
	FREE(top_errs);
	FREE(bottom_errs);
	FREE(symbols);
	FREE(channels);
      }
      /* move to the next data */
      if (plot_type & PLOT_ALL)
	cycle_data=cycle_data->next;
      else if (plot_type & PLOT_SEPARATE)
	cycle_data=cycle_data->separate_next;
    }
  }

  /* now check that we've actually output a plot */
  if (made_plot==0){
    /* we haven't actually made a plot, so if we've just output a file we should 
       delete it because it doesn't contain any information */
    if (output_medium(PGPLOT_device)==PGPLOT_MEDIUM_PAPER){
      if (strlen(opened_pgplot_device_name)!=0){
	cpgclos();
	opened_pgplot_device_name[0]='\0';
      }
      strip_pgplot_device_specifier(PGPLOT_device,tmpfile);
      remove(tmpfile);
    }
  }
  
  return(0);
}

int output_medium(char *PGPLOT_device){
  /* this routine takes a PGPLOT device specifier and returns whether it thinks
     it is a screen device or a paper device */
  int medium=PGPLOT_MEDIUM_SCREEN; /* we assume by default that it is a screen medium */

  if (strncmp(PGPLOT_device+strlen(PGPLOT_device)-2,"ps",2)==0)
    /* this is a paper output */
    medium=PGPLOT_MEDIUM_PAPER;

  return(medium);
}

void which_panel(int baseline,int if_no,int *xpanel,int *ypanel){
  /* this routine takes a baseline number and figures out which panel
     in the peak summary page it should be plotted in */
  int ant1,ant2,cx,cy,cif;
  /* the baselines are plotted like this: 
     12 13 14 15 16
     23 24 25 26 34  IF1
     35 36 45 46 56
     12 13 14 15 16
     23 24 25 26 34  IF2
     35 36 45 46 56
     xpanel=0, ypanel=0 is 35 IF2 (ie. bottom left)
     xpanel=4, ypanel=5 is 16 IF1 (ie. top right)
  */
  cx=0; cy=PANEL_Y_ROWS-1;
  for (cif=PANEL_MIN_IF;cif<=PANEL_MAX_IF;cif++){
    for (ant1=PANEL_MIN_ANT;ant1<PANEL_MAX_ANT;ant1++){
      for (ant2=ant1+1;ant2<=PANEL_MAX_ANT;ant2++){
	/* is this the baseline and IF we're after? */
	if ((baseline==(256*ant1+ant2))&&
	    (if_no==cif)){
	  /* yes it is */
	  *xpanel=cx;
	  *ypanel=cy;
	  return; /* no need to do any more */
	}
	cx++;
	if (cx==PANEL_X_COLS){
	  cx=0; cy--;
	}
      }
    }
  }
}

int add_correlation(correlation_data **correlation_list,float **vis,float **wgt,int baseline,
		    float ut,float u,float v,float w,int flag,int bin,int if_no,int sourceno,
		    int if_nfreq,int if_nstok,int longest_baseline,int scannum){
  /* this routine adds a correlation point to the correlation_data linked list at the
     end of the list */
  int i,nodata=0,wgt_size;
  correlation_data *cycle_data=NULL,*new_correlation=NULL;
  struct tm obs_time;
  char datobs[BUFSIZE];

  /* we don't add the correlation if the data is flagged as bad */
  if (flag!=0)
    return(-1);

  /* this loop gets us to the last entry in the list */
  if (*correlation_list==NULL){
    /* there is no data currently */
    nodata=1;
    cycle_data=NULL;
  } else 
    /* move cycle_data to the end of the list */
    for (cycle_data=*correlation_list;cycle_data->next!=NULL;cycle_data=cycle_data->next);
  
  /* allocate memory for our new entry */
  MALLOC(new_correlation,sizeof(correlation_data));
  wgt_size=if_nstok*if_nfreq;
  MALLOC(new_correlation->vis,2*wgt_size*sizeof(float));
  MALLOC(new_correlation->wgt,wgt_size*sizeof(float));

  /* initialise some pointers */
  new_correlation->spectrum_info.peak_channels=NULL;
  new_correlation->spectrum_info.peak_verified=NULL;
  new_correlation->spectrum_info.peak_number=NULL;
  new_correlation->spectrum_info.peak_present_both_pols=NULL;
  new_correlation->individual_pols=NULL;

  /* now assign the values */
  /* the baseline number, and whether this is the longest baseline 
     in the array */
  new_correlation->baseline=baseline;
  if (baseline==longest_baseline)
    new_correlation->islongest_baseline=1;
  else
    new_correlation->islongest_baseline=0;

  /* the bandwidth of the observation */
  for (i=0;i<if_.n_if;i++){
    if (if_.if_num[i]==if_no){
      new_correlation->if_bandwidth=doubles_.if_bw[i];
      break;
    }
  }

  /* the time of the observation */
  new_correlation->ut=ut;
  obs_time.tm_sec=(int)ut; /* coz ut is given in seconds since midnight */
  obs_time.tm_min=obs_time.tm_hour=0;
  strncpy(datobs,names_.datobs,12);
  datobs[12]='\0';
  strip_trailing_spaces(datobs);
  sscanf(datobs,"%d-%2d-%2d",&(obs_time.tm_year),&(obs_time.tm_mon),&(obs_time.tm_mday));
  if (obs_time.tm_year>100){
    /* four digit year, need to subtract 1900 to fit into C's time model */
    obs_time.tm_year-=1900;
  } /* else the year is a two digit year, which is already the number of years since 1900 */
  obs_time.tm_mon-=1; /* months since January */
  new_correlation->observation_time=timegm(&obs_time); /* convert to simple time representation */
  
  /* the uvw plane location */
  new_correlation->u=u;
  new_correlation->v=v;
  new_correlation->w=w;

  /* is the data alright? */
  new_correlation->flag=flag;
  new_correlation->bin=bin;
/*   new_correlation->polarisation_non_detection=0; */

  /* the IF number */
  new_correlation->if_no=if_no;

  /* details of the source being observed */
  new_correlation->sourceno=sourceno;
  for (i=0;i<su_.n_su;i++){
    if (sourceno==su_.su_num[i]){
      strncpy(new_correlation->sourcename,names_.su_name+i*16,16);
      new_correlation->sourcename[15]='\0';
      strip_trailing_spaces(new_correlation->sourcename);
      break;
    }
  }

  /* the names of the polarisations that were observed */
  new_correlation->stokes_name=NULL;
  MALLOC(new_correlation->stokes_name,if_nstok*sizeof(char*));
  for (i=0;i<if_nstok;i++){
    MALLOC(new_correlation->stokes_name[i],BUFSIZE*sizeof(char));
    strncpy(new_correlation->stokes_name[i],names_.if_cstok+(if_no-1)*if_nstok*2+i*2,2);
    new_correlation->stokes_name[i][2]='\0';
  }

  /* the number of frequency channels */
  new_correlation->if_nfreq=if_nfreq;

  /* the number of Stokes parameters observed */
  new_correlation->if_nstok=if_nstok;

  /* the number of the scan since the observation started */
  new_correlation->scannum=scannum;

  /* the data and the weight assigned to it */
  for (i=0;i<wgt_size;i++){
    new_correlation->wgt[i]=(*wgt)[i];
    new_correlation->vis[2*i]=(*vis)[2*i];
    new_correlation->vis[2*i+1]=(*vis)[2*i+1];
  }
  new_correlation->pols_combined=NULL;

  /* now attach it to the the end of the list */
  new_correlation->next=NULL;
  new_correlation->prev=cycle_data;
  new_correlation->separate_next=NULL;
  if (nodata==1)
    *correlation_list=new_correlation;
  else
    cycle_data->next=new_correlation;

  /* connect the baselines */
  for (cycle_data=new_correlation->prev;cycle_data!=NULL;
       cycle_data=cycle_data->prev){
    if ((cycle_data->baseline==new_correlation->baseline)&&
	(cycle_data->if_no==new_correlation->if_no)){
      cycle_data->samebaseline_next=new_correlation;
      new_correlation->samebaseline_prev=cycle_data;
      break;
    }
  }

  return(0);
}

int clear_correlations(correlation_data **all_correlations){
  /* this routine clears all the memory in the correlation list */
  correlation_data *free_correlation=NULL;

  while(*all_correlations!=NULL){
    free_correlation=*all_correlations;
    *all_correlations=(*all_correlations)->next;
    remove_correlation(&free_correlation);
  }

  return(0);
}

int remove_correlation(correlation_data **correlation_to_remove){
  /* this frees the memory used by a correlation element */
  int i;
  FREE((*correlation_to_remove)->vis);
  FREE((*correlation_to_remove)->wgt);
  for (i=0;i<(*correlation_to_remove)->if_nstok;i++){
    FREE((*correlation_to_remove)->stokes_name[i]);
    if ((*correlation_to_remove)->individual_pols!=NULL)
      FREE((*correlation_to_remove)->individual_pols[i]);
  }
  FREE((*correlation_to_remove)->stokes_name);
  FREE((*correlation_to_remove)->individual_pols);
  FREE((*correlation_to_remove)->pols_combined);
  FREE((*correlation_to_remove)->spectrum_info.peak_channels);
  FREE((*correlation_to_remove)->spectrum_info.peak_verified);
  FREE((*correlation_to_remove)->spectrum_info.peak_number);
  FREE(*correlation_to_remove);

  return(0);
}

int separate_lists(correlation_data **all_correlations,
		   correlation_data ***separate_correlations){
  /* this routine separates the all_correlations list of data into a number 
     of lists, one for each source name */
  int n_lists=0,need_new_list,i;
  correlation_data *cycle_data=NULL,*cycle_separate=NULL;
  
  for (cycle_data=*all_correlations;cycle_data!=NULL;cycle_data=cycle_data->next){
    /* check if a list exists already for this source name */
    need_new_list=1;
    for (i=0;i<n_lists;i++){
      if (strcmp(((*separate_correlations)[i])->sourcename,cycle_data->sourcename)==0){
	/* we have a source name match */
	need_new_list=0;
	/* add this data to the end of this list */
	for (cycle_separate=(*separate_correlations)[i];cycle_separate->separate_next!=NULL;
	     cycle_separate=cycle_separate->separate_next);
	cycle_separate->separate_next=cycle_data;
	cycle_data->separate_next=NULL;
	break;
      }
    }
    if (need_new_list==1){
      /* need to make more space for a new list */
      n_lists++;
      *separate_correlations=realloc(*separate_correlations,n_lists*sizeof(correlation_data*));
      (*separate_correlations)[n_lists-1]=cycle_data;
      ((*separate_correlations)[n_lists-1])->separate_next=NULL;
    }
  }

  return(n_lists);
}


void palett(int type,float contra,float bright){
  /* this routine sets up a colour mapping for the PGPLOT device and was 
     shamelessly ripped off from pgdemo4 */
  static float gl[]={0.0, 1.0};
  static float gr[]={0.0, 1.0};
  static float gg[]={0.0, 1.0};
  static float gb[]={0.0, 1.0};

  static float rl[]={-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7};
  static float rr[]={ 0.0, 0.0, 0.0 , 0.0 , 0.6 , 1.0 , 1.0 , 1.0, 1.0};
  static float rg[]={ 0.0, 0.0, 0.0 , 1.0 , 1.0 , 1.0 , 0.6 , 0.0, 1.0};
  static float rb[]={ 0.0, 0.3, 0.8 , 1.0 , 0.3 , 0.0 , 0.0 , 0.0, 1.0};

  static float hl[]={0.0, 0.2, 0.4, 0.6, 1.0};
  static float hr[]={0.0, 0.5, 1.0, 1.0, 1.0};
  static float hg[]={0.0, 0.0, 0.5, 1.0, 1.0};
  static float hb[]={0.0, 0.0, 0.0, 0.3, 1.0};

  static float wl[]={0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0};
  static float wr[]={0.0, 1.0, 0.0, 0.0, 0.3, 0.8 , 0.3 , 1.0 , 1.0 , 1.0};
  static float wg[]={0.0, 0.5, 0.4, 1.0, 0.0, 0.0 , 0.2 , 0.7 , 1.0 , 1.0};
  static float wb[]={0.0, 0.0, 0.0, 0.0, 0.4, 1.0 , 0.0 , 0.0 , 0.95, 1.0};

  static float al[]={0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5,
		     0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0};
  static float ar[]={0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 
		     0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  static float ag[]={0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8,
		     0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0};
  static float ab[]={0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9,
		     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  if (type==PGPLOT_CT_GREYSCALE) /* greyscale */
    cpgctab(gl,gr,gg,gb,2,contra,bright);
  else if (type==PGPLOT_CT_RAINBOW) /* rainbow */
    cpgctab(rl,rr,rg,rb,9,contra,bright);
  else if (type==PGPLOT_CT_HEAT) /* heat */
    cpgctab(hl,hr,hg,hb,5,contra,bright);
  else if (type==PGPLOT_CT_WEIRDIRAF) /* weird IRAF */
    cpgctab(wl,wr,wg,wb,10,contra,bright);
  else if (type==PGPLOT_CT_AIPS) /* AIPS */
    cpgctab(al,ar,ag,ab,20,contra,bright);
}

void setvp(void){
  /* another routine shamelessly lifted from pgdemo4, this time to set
     the viewport */
  float d,vpx1,vpx2,vpy1,vpy2;

  cpgsvp(0.0,1.0,0.0,1.0);
  cpgqvp(1,&vpx1,&vpx2,&vpy1,&vpy2);
  if ((vpx2-vpx1)<(vpy2-vpy1))
    d=(vpx2-vpx1)/40.0;
  else
    d=(vpy2-vpy1)/40.0;
  vpx1+=5.0*d;
  vpx2-=2.0*d;
  vpy1+=8.0*d;
  vpy2-=2.0*d;
  cpgvsiz(vpx1,vpx2,vpy1,vpy2);
}

int change_to_lagspace(correlation_data **all_correlations,
		       correlation_data **all_transforms,
		       int cycle_type){
  /* this routine takes the all_correlations list and make a new list which 
     will be the fourier transform of the all_correlations visibility */
  correlation_data *new_transform=NULL,*cycle_data=NULL;
  correlation_data *cycle_transforms=NULL,*second_cycle=NULL;
  fftw_complex *in=NULL,*out=NULL;
  fftw_plan p;
  int i,j;

  cycle_data=*all_correlations;
  if (cycle_data==NULL)
    return(0);

  /* prepare the FFT */
  if (cycle_data!=NULL){
    in=fftw_malloc(cycle_data->if_nfreq*sizeof(fftw_complex));
    out=fftw_malloc(cycle_data->if_nfreq*sizeof(fftw_complex));
    p=fftw_plan_dft_1d(cycle_data->if_nfreq,in,out,FFTW_BACKWARD,FFTW_ESTIMATE);
  }
  
  /* now go through the list */
  while(cycle_data!=NULL){
    /* prepare a new transform */
    MALLOC(new_transform,sizeof(correlation_data));

    /* copy the basic information over */
    copy_correlation_info(&cycle_data,&new_transform);
    /* copy the values to the in array for each polarisation */
    for (i=0;i<cycle_data->if_nstok;i++){
      for (j=0;j<cycle_data->if_nfreq;j++){
	in[j]=cycle_data->vis[j*2*cycle_data->if_nstok+i*2+0]+
	  cycle_data->vis[j*2*cycle_data->if_nstok+i*2+1]*I;
      }
      /* now do the FFT */
      fftw_execute(p);
      /* now copy the out array to the new structure */
      for (j=0;j<cycle_data->if_nfreq;j++){
	new_transform->vis[j*2*cycle_data->if_nstok+i*2+0]=crealf(out[j])/sqrt(cycle_data->if_nfreq);
	new_transform->vis[j*2*cycle_data->if_nstok+i*2+1]=cimagf(out[j])/sqrt(cycle_data->if_nfreq);
      }
    }
    /* add the new transform to the list */
    new_transform->next=NULL;
    new_transform->prev=NULL;
    new_transform->samebaseline_next=NULL;
    new_transform->samebaseline_prev=NULL;
    new_transform->separate_next=NULL;
    if (*all_transforms!=NULL){
      if (cycle_type & PLOT_ALL){
	for (cycle_transforms=*all_transforms;cycle_transforms->next!=NULL;
	     cycle_transforms=cycle_transforms->next);
	cycle_transforms->next=new_transform;
	new_transform->prev=cycle_transforms;
      } else if (cycle_type & PLOT_SEPARATE) {
	for (cycle_transforms=*all_transforms;cycle_transforms->separate_next!=NULL;
	     cycle_transforms=cycle_transforms->separate_next);
	cycle_transforms->separate_next=new_transform;
	new_transform->prev=cycle_transforms;
      }
    } else {
      *all_transforms=new_transform;
    }
    /* connect the baselines */
    for (second_cycle=new_transform->prev;second_cycle!=NULL;
	 second_cycle=second_cycle->prev){
      if ((second_cycle->baseline==new_transform->baseline)&&
	  (second_cycle->if_no==new_transform->if_no)){
	second_cycle->samebaseline_next=new_transform;
	new_transform->samebaseline_prev=second_cycle;
	break;
      }
    }
    if (cycle_type & PLOT_ALL)
      cycle_data=cycle_data->next;
    else if (cycle_type & PLOT_SEPARATE)
      cycle_data=cycle_data->separate_next;
  }

  /* deallocate the plan and free our memory */
  fftw_destroy_plan(p);
  fftw_free(in);
  in=NULL;
  fftw_free(out);
  out=NULL;

  return(0);
}

#define COPY_INFO_NUMBER(a,b,c) (*a)->c=(*b)->c
#define COPY_INFO_STRING(a,b,c) strcpy((*a)->c,(*b)->c)

void copy_correlation_info(correlation_data **original,
			   correlation_data **destination){
  /* this routine copies over the basic information (ie. everything but the vis array
     and the list pointers) from one correlation_data structure to another */
  int i,wgtsize;

  /* copy all the simple data */
  COPY_INFO_NUMBER(destination,original,baseline);
  COPY_INFO_NUMBER(destination,original,islongest_baseline);
  COPY_INFO_NUMBER(destination,original,ut);
  COPY_INFO_NUMBER(destination,original,observation_time);
  COPY_INFO_NUMBER(destination,original,u);
  COPY_INFO_NUMBER(destination,original,v);
  COPY_INFO_NUMBER(destination,original,w);
  COPY_INFO_NUMBER(destination,original,flag);
  COPY_INFO_NUMBER(destination,original,bin);
  COPY_INFO_NUMBER(destination,original,if_no);
  COPY_INFO_NUMBER(destination,original,sourceno);
  COPY_INFO_NUMBER(destination,original,if_nfreq);
  COPY_INFO_NUMBER(destination,original,if_nstok);
  COPY_INFO_NUMBER(destination,original,if_bandwidth);
  COPY_INFO_NUMBER(destination,original,scannum);
/*   COPY_INFO_NUMBER(destination,original,polarisation_non_detection); */
  COPY_INFO_STRING(destination,original,sourcename);

  /* initialise some pointers */
  (*destination)->spectrum_info.peak_channels=NULL;
  (*destination)->spectrum_info.peak_verified=NULL;
  (*destination)->spectrum_info.peak_number=NULL;
  (*destination)->spectrum_info.peak_present_both_pols=NULL;
  (*destination)->individual_pols=NULL;

  /* copy the Stokes name data */
  MALLOC((*destination)->stokes_name,(*destination)->if_nstok*sizeof(char*));
  for (i=0;i<(*destination)->if_nstok;i++){
    MALLOC((*destination)->stokes_name[i],BUFSIZE*sizeof(char));
    COPY_INFO_STRING(destination,original,stokes_name[i]);
  }

  /* copy the weighting data */
  wgtsize=((*original)->if_nfreq) * ((*original)->if_nstok);
  (*destination)->wgt=NULL;
  (*destination)->vis=NULL;
  MALLOC((*destination)->wgt,wgtsize*sizeof(float));
  for (i=0;i<wgtsize;i++)
    COPY_INFO_NUMBER(destination,original,wgt[i]);
  /* allocate the vis array */
  MALLOC((*destination)->vis,2*wgtsize*sizeof(float));

}

int combine_polarisations(correlation_data **all_correlations,int cycle_type,int stokes_to_make){
  /* this routine goes through the list all_correlations and scalar adds the two
     polarisations together */
  correlation_data *cycle_data=NULL;
  int i,j;
  float average_amplitude,rms_noise;

  cycle_data=*all_correlations;
  if (cycle_data==NULL)
    return(0);

  while(cycle_data!=NULL){
    /* which Stokes parameter are we making */
    if (stokes_to_make==STOKES_I){
      /* Stokes I is just (XX + YY)/2 */
      /* allocate the necessary memory */
      CALLOC(cycle_data->pols_combined,cycle_data->if_nfreq,sizeof(float));
      MALLOC(cycle_data->individual_pols,cycle_data->if_nstok*sizeof(float*));
      /* go through the Stokes parameters and search for the ones we want */
      for (i=0;i<cycle_data->if_nstok;i++){
	if ((strcmp(cycle_data->stokes_name[i],"XX")==0)||
	    (strcmp(cycle_data->stokes_name[i],"YY")==0)){
	  MALLOC(cycle_data->individual_pols[i],cycle_data->if_nfreq*sizeof(float));
	  for (j=0;j<cycle_data->if_nfreq;j++){
	    /* alright, we should note that if you combine the pols in this way after
	       doing the Fourier transforms, then actually you're not dealing with 
	       XX and YY anymore, so can't make Stokes I with (XX + YY)/2. Instead
	       what we are trying to do is reduce the amount of noise in the spectrum
	       by adding the two individual "polarisations" together, so we use the
	       square root of the squares */
	    cycle_data->pols_combined[j]+=powf(cycle_data->vis[j*2*cycle_data->if_nstok+i*2],2); /* we only take the real bit */
	    cycle_data->individual_pols[i][j]=fabsf(cycle_data->vis[j*2*cycle_data->if_nstok+i*2]); /* we make it positive */
	  }
	} else 
	  cycle_data->individual_pols[i]=NULL;
      }
      for (j=0;j<cycle_data->if_nfreq;j++)
	cycle_data->pols_combined[j]=sqrtf(cycle_data->pols_combined[j]);
      strcpy(cycle_data->combined_pol_name,"I");
    } else if (stokes_to_make==STOKES_Q){
      ERROR("making Stokes Q is currently not supported\n");
    } else if (stokes_to_make==STOKES_U){
      ERROR("making Stokes U is currently not supported\n");
    } else if (stokes_to_make==STOKES_V){
      ERROR("making Stokes V is currently not supported\n");
    }

    /* now we compute the RMS noise and average amplitude of the spectrum we just made */
    rms_noise=standard_deviation(&(cycle_data->pols_combined),
				 cycle_data->if_nfreq,&average_amplitude);
    cycle_data->spectrum_info.average_amplitude=average_amplitude;
    cycle_data->spectrum_info.rms_noise=rms_noise;
    cycle_data->spectrum_info.npeaks_found=0;
    cycle_data->spectrum_info.peak_channels=NULL;
    cycle_data->spectrum_info.peak_verified=NULL;
    cycle_data->spectrum_info.peak_number=NULL;

    if (cycle_type & PLOT_ALL)
      cycle_data=cycle_data->next;
    else if (cycle_type & PLOT_SEPARATE)
      cycle_data=cycle_data->separate_next;
  }

  return(0);
}

void shift_data(float *data,int ndata,int nshift){
  float *newdata=NULL;
  int i,j;
  
  CALLOC(newdata,ndata,sizeof(float));
  for (i=0;i<ndata;i++){
    j=i+nshift;
    if (j>=ndata)
      j-=ndata;
    newdata[j]=data[i];
  }
  for (i=0;i<ndata;i++){
    data[i]=newdata[i];
  }
  FREE(newdata);
}

void findpeaks(float *data,int ndata,float slopethreshold,float ampthreshold,
	       float smoothwidth,int peakgroup,int *npeaks,int **peakchannel,
	       float **peakheight,float **peakwidth){
  /* 
     Function to locate the positive peaks in a noisy x-y data set.
     Detects peaks by looking for downward zero-crossings in the
     first derivative that exceed slopethreshold. Returns npeaks
     as the number of peaks it has found, along with the channel it
     was found in, its height and width. slopethreshold, ampthreshold
     and smoothwidth control sensitivity. Higher values will neglect
     smaller features. peakgroup is the number of points around the
     "top part" of the peak.
   */

  float tmpf,*smoothed_data=NULL,n,AmpTest,*derivdata=NULL;
  float *yy=NULL,PeakY,*values=NULL,*shiftdata=NULL;
  int j,groupindex,k,PeakX,i,*indices=NULL,nindices,shift;
  int *po=NULL;

  tmpf=roundf(smoothwidth);
  smoothwidth=tmpf;
  
  CALLOC(po,ndata,sizeof(int));
  CALLOC(shiftdata,ndata,sizeof(float));
  for (shift=0;shift<ndata;shift++){
    /* make a copy of the data first */
    for (i=0;i<ndata;i++)
      shiftdata[i]=data[i];
    /* shift by the shift number */
    shift_data(shiftdata,ndata,shift);
    
    /* take the derivative */
    deriv(shiftdata,ndata,&derivdata);

    fastsmooth(derivdata,ndata,(int)smoothwidth,&smoothed_data);
    n=floorf((float)peakgroup/2.0);
    *npeaks=0;
    AmpTest=ampthreshold;
    MALLOC(yy,peakgroup*sizeof(float));
    for (j=0;j<(ndata-(int)smoothwidth);j++){
      if (((smoothed_data[j]>0)-(smoothed_data[j]<0))>
	  ((smoothed_data[j+1]>0)-(smoothed_data[j+1]<0))){
	/* detects zero crossing */
	if ((smoothed_data[j]-smoothed_data[j+1])>slopethreshold*data[j]){
	  /* if slope of derivative is larger than slopethreshold */
	  if (data[j]>AmpTest){
	    /* if height of peak is larger than ampthreshold */
	    for (k=0;k<peakgroup;k++){
	      groupindex=j+k-n;
	      if (groupindex<0)
		groupindex=0;
	      if (groupindex>=ndata)
		groupindex=ndata;
	      yy[k]=shiftdata[groupindex];
	    }
	    PeakY=yy[0];
	    for (i=1;i<peakgroup;i++){
	      if (yy[i]>PeakY)
		PeakY=yy[i];
	    }
	    val2ind(data,ndata,PeakY,&indices,&values,&nindices);
	    PeakX=indices[0];
	    po[PeakX]++;
	  }
	}
      }
    }
    FREE(yy);
  }
  FREE(shiftdata);

  for (i=0;i<ndata;i++){
    if (po[i]>1){
      (*npeaks)++;
      REALLOC(*peakchannel,*npeaks*sizeof(int));
      (*peakchannel)[*npeaks-1]=i;
      REALLOC(*peakheight,*npeaks*sizeof(float));
      (*peakheight)[*npeaks-1]=data[i];
      REALLOC(*peakwidth,*npeaks*sizeof(float));
      (*peakwidth)[*npeaks-1]=0.0;
    }
  }
  FREE(po);
}

void fastsmooth(float *data,int ndata,int smoothwidth,float **smoothed_data){
  float SumPoints,*s=NULL,*SmoothData=NULL;
  int i,halfw,k;

  for (i=0,SumPoints=0.0;i<smoothwidth;i++)
    SumPoints+=data[i];
  MALLOC(s,ndata*sizeof(float));
  halfw=(int)(roundf((float)smoothwidth/2.0));
  for (k=0;k<ndata-smoothwidth;k++){
    s[k+halfw-1]=SumPoints;
    SumPoints-=data[k];
    SumPoints+=data[k+smoothwidth];
  }
  for (i=0;i<ndata;i++){
    s[i]=s[i]/(float)smoothwidth;
  }

  for (i=0,SumPoints=0.0;i<smoothwidth;i++)
    SumPoints+=s[i];
  SmoothData=calloc(ndata,sizeof(float));

  for (k=0;k<ndata-smoothwidth;k++){
    SmoothData[k+halfw-1]=SumPoints;
    SumPoints-=s[k];
    SumPoints+=s[k+smoothwidth];
  }
  for (i=0;i<ndata;i++){
    SmoothData[i]=SmoothData[i]/(float)smoothwidth;
  }

  FREE(s);
  *smoothed_data=SmoothData;
}

void deriv(float *data,int ndata,float **deriv_data){
  /*
    First derivative of vector using 2-point central difference.
    Example: deriv([1 1 1 2 3 4]) yields ([0 0 0.5 1 1 1])
   */

  float *d=NULL;
  int j;

  d=calloc(ndata,sizeof(float));
  d[0]=data[1]-data[0];
  d[ndata-1]=data[ndata-1]-data[ndata-2];
  for (j=1;j<ndata-1;j++){
    d[j]=(data[j+1]-data[j-1])/2.0;
  }
  *deriv_data=d;
}

void val2ind(float *data,int ndata,float val,int **indices,float **values,
	     int *nindices){
  /* 
     Returns the index of the element of vector x that is closest to
     val.
   */
  float dif,mindif;
  int i;

  *nindices=0;

  mindif=fabsf(data[0]-val);
  for (i=1;i<ndata;i++){
    dif=fabsf(data[i]-val);
    if (dif<mindif){
      mindif=dif;
    }
  }

  for (i=0;i<ndata;i++){
    dif=fabsf(data[i]-val);
    if (dif==mindif){
      (*nindices)++;
      *indices=realloc(*indices,*nindices*sizeof(int));
      *values=realloc(*values,*nindices*sizeof(float));
      (*indices)[*nindices-1]=i;
      (*values)[*nindices-1]=data[i];
    }
  }

}

int peak_finder(correlation_data **all_transforms,correlation_data ***peak_lists){
  /* this routine takes the all_correlations list, sorts it into time order, 
     then looks for peaks in the lag spectra, and returns each peak as a list of
     2N+1 spectra, where spectra 1 -> N are before the peak, spectrum N+1 is the
     spectrum in which the peak is found, and spectra N+1 -> 2N+1 are the spectra
     after the peak */
  correlation_data *cycle_data=NULL,*second_cycle=NULL;
  correlation_data **nearby_data_before=NULL,**nearby_data_after=NULL;
  float *rms_before=NULL,*rms_after=NULL,average_rms_before,average_rms_after;
  float *average_before=NULL,*average_after=NULL;
  float average_rms,min_peak_level;
  float slopethreshold,ampthreshold,smoothwidth,*peakheight=NULL,*peakwidth=NULL;
  int i,nbefore,nafter,nfitpoints,j,widthpoints;
  int npeaks,*peakchannel=NULL;
  double c0,c1,cov00,cov01,cov11,chisq;
  double *xpoints=NULL,*ypoints=NULL,*weights=NULL;
  double average_estimate,average_estimate_error;

  if (*all_transforms==NULL)
    return(0);

  /* first thing to do is make sure that the list is in time order */
  time_sort_data(all_transforms);

  /* we're going to keep a list of the N spectra before and after each 
     spectrum in the list */
  MALLOC(nearby_data_before,arguments.number_of_nearby_spectra*sizeof(correlation_data*));
  MALLOC(nearby_data_after,arguments.number_of_nearby_spectra*sizeof(correlation_data*));
  /* and the RMS noise values for those spectra */
  MALLOC(rms_before,arguments.number_of_nearby_spectra*sizeof(float));
  MALLOC(rms_after,arguments.number_of_nearby_spectra*sizeof(float));
  /* and the average baseline value for those spectra */
  MALLOC(average_before,arguments.number_of_nearby_spectra*sizeof(float));
  MALLOC(average_after,arguments.number_of_nearby_spectra*sizeof(float));

  /* now traverse the list */
  cycle_data=*all_transforms;
  while(cycle_data!=NULL){
    /* clear the nearby data array */
    for (i=0;i<arguments.number_of_nearby_spectra;i++){
      nearby_data_before[i]=NULL;
      nearby_data_after[i]=NULL;
    }
    nbefore=0;
    nafter=0;
    /* find the locations of the nearby spectra */
    /* first find the earlier spectra, remembering we need to keep to the
       same baseline */
    for (i=0,second_cycle=cycle_data->samebaseline_prev;
	 ((i<arguments.number_of_nearby_spectra)&&(second_cycle!=NULL));
	 i++,second_cycle=second_cycle->samebaseline_prev){
      nearby_data_before[i]=second_cycle;
    }
    nbefore=i;
    /* then find the later spectra */
    for (i=0,second_cycle=cycle_data->samebaseline_next;
	 ((i<arguments.number_of_nearby_spectra)&&(second_cycle!=NULL));
	 i++,second_cycle=second_cycle->samebaseline_next){
      nearby_data_after[i]=second_cycle;
    }
    nafter=i;

    /* now get the RMS noise figure for the spectra before and after 
       this spectrum */
    for (i=0;i<nbefore;i++){
      /* get the standard deviation for the spectra, which is also the RMS
         noise level */
      rms_before[i]=(nearby_data_before[i])->spectrum_info.rms_noise;
      average_before[i]=(nearby_data_before[i])->spectrum_info.average_amplitude;
    }
    for (i=0;i<nafter;i++){
      rms_after[i]=(nearby_data_after[i])->spectrum_info.rms_noise;
      average_after[i]=(nearby_data_after[i])->spectrum_info.average_amplitude;
    }
    /* get the average values of these RMS */
    average_rms_before=meanvalue(&rms_before,nbefore);
    average_rms_after=meanvalue(&rms_after,nafter);
    /* and the average of those two values */
    average_rms=(average_rms_before+average_rms_after)/2.0;
    cycle_data->spectrum_info.estimated_rms_noise=average_rms;

    /* we need to know the average level of the spectrum as well, but of course
       we can't use the measured value because a large peak may skew the average.
       so we compute what the average might be by taking a linear fit through
       the nearby average values. we'll use the gnu scientific library fitting
       functions for this */
    /* prepare the arrays we need to pass to the fitting routine */
    nfitpoints=nbefore+nafter;
    MALLOC(xpoints,nfitpoints*sizeof(double));
    MALLOC(ypoints,nfitpoints*sizeof(double));
    MALLOC(weights,nfitpoints*sizeof(double));
    
    for (i=0,j=0;i<nbefore;i++,j++){
      xpoints[j]=(nearby_data_before[i])->observation_time-cycle_data->observation_time;
      ypoints[j]=average_before[i];
      /* assign the weights, as the reciprocal of the square of the standard deviation
	 of the point */
      weights[j]=1.0/pow(rms_before[i],2);
    }
    for (i=0;i<nafter;i++,j++){
      xpoints[j]=(nearby_data_after[i])->observation_time-cycle_data->observation_time;
      ypoints[j]=average_after[i];
      weights[j]=1.0/pow(rms_after[i],2);
    }

    /* now make the fit */
    gsl_fit_wlinear(xpoints,1,weights,1,ypoints,1,nfitpoints,
		    &c0,&c1,&cov00,&cov01,&cov11,&chisq);

    /* use the fit to calculate the likely average for this spectrum */
    gsl_fit_linear_est(0,c0,c1,cov00,cov01,cov11,
		       &average_estimate,&average_estimate_error);
    cycle_data->spectrum_info.estimated_average_amplitude=average_estimate;
    cycle_data->spectrum_info.estimated_average_amplitude_error=average_estimate_error;
    cycle_data->spectrum_info.fit_c0=c0;
    cycle_data->spectrum_info.fit_c1=c1;

    /* now we pass the spectrum to the findpeaks routine */
    widthpoints=10;
    slopethreshold=0.5*powf(widthpoints,-2);
    ampthreshold=cycle_data->spectrum_info.estimated_average_amplitude*0.1;
    smoothwidth=3;
    findpeaks(cycle_data->pols_combined,cycle_data->if_nfreq,
	       slopethreshold,ampthreshold,smoothwidth,5,&npeaks,
	       &peakchannel,&peakheight,&peakwidth);
    /* we now have a list of peaks in the spectrum */
    min_peak_level=average_estimate+arguments.minimum_peak_sigma_factor
      *average_rms;
    for (i=0;i<npeaks;i++){
      /* check the peak satisfies the following criteria */
      /* is the peak height more than N times the RMS noise level? */
      if (peakheight[i]<min_peak_level){
	/* no! */
	continue;
      }
      /* are we looking at the longest baseline? */
      if (cycle_data->islongest_baseline==1){
	/* yes, so we ignore peaks in channel 0. it is probably just
	   a continuum source at the phase centre, and even if it is a
	   transient, we can't tell these sources apart from continuum
	   sources, so we just don't bother */
	if (peakchannel[i]==0){
	  continue;
	}
      }
      /* this is a peak worthy of further examination */
      cycle_data->spectrum_info.npeaks_found++;
      REALLOC(cycle_data->spectrum_info.peak_channels,
	      cycle_data->spectrum_info.npeaks_found*sizeof(int));
      REALLOC(cycle_data->spectrum_info.peak_verified,
	      cycle_data->spectrum_info.npeaks_found*sizeof(int));
      REALLOC(cycle_data->spectrum_info.peak_number,
	      cycle_data->spectrum_info.npeaks_found*sizeof(int));
      REALLOC(cycle_data->spectrum_info.peak_present_both_pols,
		cycle_data->spectrum_info.npeaks_found*sizeof(int));
      cycle_data->spectrum_info.peak_channels
	[cycle_data->spectrum_info.npeaks_found-1]=peakchannel[i];
      cycle_data->spectrum_info.peak_verified
	[cycle_data->spectrum_info.npeaks_found-1]=PEAK_VERIFIED_NO;
      cycle_data->spectrum_info.peak_number
	[cycle_data->spectrum_info.npeaks_found-1]=-1;
      cycle_data->spectrum_info.peak_present_both_pols
	[cycle_data->spectrum_info.npeaks_found-1]=1;
    }

/*     /\* now we look at the current spectra, and we find any peak that has an */
/*        amplitude of N times this average RMS noise level on top of the average */
/*        level that we've estimated for this spectrum *\/ */
/*     min_peak_level=average_estimate+arguments.minimum_peak_sigma_factor*average_rms; */
/*     if (cycle_data->islongest_baseline==1){ */
/*       /\* if there is a peak in channel 0 on the longest baseline, then we */
/*          are going to assume that a continuum source is at the phase centre, */
/*          so we will not mark such peaks *\/ */
/*       start_channel=1; */
/*     } else { */
/*       /\* we don't ignore peaks in channel 0 for baselines other than the */
/*          longest one, since we were missing a lot of potential peaks that */
/*          were scaling into the first lag channel *\/ */
/*       start_channel=0; */
/*     } */
/*     for (i=start_channel;i<cycle_data->if_nfreq;i++){ */
/*       if (cycle_data->pols_combined[i]>=min_peak_level){ */
/* 	/\* we've found a peak that is big enough to warrant further examination *\/ */
/* 	cycle_data->spectrum_info.npeaks_found++; */
/* 	REALLOC(cycle_data->spectrum_info.peak_channels, */
/* 		cycle_data->spectrum_info.npeaks_found*sizeof(int)); */
/* 	REALLOC(cycle_data->spectrum_info.peak_verified, */
/* 		cycle_data->spectrum_info.npeaks_found*sizeof(int)); */
/* 	REALLOC(cycle_data->spectrum_info.peak_number, */
/* 		cycle_data->spectrum_info.npeaks_found*sizeof(int)); */
/* 	REALLOC(cycle_data->spectrum_info.peak_present_both_pols, */
/* 		cycle_data->spectrum_info.npeaks_found*sizeof(int)); */
/* 	cycle_data->spectrum_info.peak_channels[cycle_data->spectrum_info.npeaks_found-1]=i; */
/* 	cycle_data->spectrum_info.peak_verified[cycle_data->spectrum_info.npeaks_found-1]= */
/* 	  PEAK_VERIFIED_NO; */
/* 	cycle_data->spectrum_info.peak_number[cycle_data->spectrum_info.npeaks_found-1]=-1; */
/* 	cycle_data->spectrum_info.peak_present_both_pols[cycle_data->spectrum_info.npeaks_found-1]=1; */
/*       } */
/*     } */
    /* free the fitting memory */
    FREE(xpoints);
    FREE(ypoints);
    FREE(weights);
    FREE(peakchannel);
    FREE(peakheight);
    FREE(peakwidth);

    /* go to the next spectrum */
    cycle_data=cycle_data->next;
  }

  /* free our memory */
  FREE(nearby_data_before);
  FREE(nearby_data_after);
  FREE(rms_before);
  FREE(rms_after);
  FREE(average_before);
  FREE(average_after);
  
  return(0);
}

int cull_peaks(correlation_data **all_transforms){
  /* this routine takes the all_transforms list after the peak finder has run on it,
     and goes through the peaks that were found, and removing those peaks that do not
     satisfy more stringent criteria */
  correlation_data *cycle_data=NULL,*second_cycle=NULL,*baseline_cycle=NULL;
  int i,j,k,found_corresponding_peak=0,found_peak_match,peak_number=0,found_baseline=0;
  int matched_baseline=0,ant1,ant2;
  float long_baseline_length,this_baseline_length,baseline_ratio;
  float long_baseline_lag,long_baseline_lag_err,this_baseline_lag,this_baseline_lag_err;
  float this_baseline_expected_lag,this_baseline_expected_lag_err;
  char timestring[BUFSIZE],tmp[BUFSIZE];

  if (*all_transforms==NULL)
    return(0);

  cycle_data=*all_transforms;
  while(cycle_data!=NULL){
    /* first we look at the peaks only in IF 1 on the longest baseline */
    if ((cycle_data->if_no==1)&&(cycle_data->islongest_baseline==1)){
      if (cycle_data->spectrum_info.npeaks_found>0){
	/* we have peaks in this spectrum */
	for (i=0;i<cycle_data->spectrum_info.npeaks_found;i++){
	  /* we don't report about peaks that are already in lag channel 0 in the
	     longest baseline! then all the calibrators show up! */
	  if (cycle_data->spectrum_info.peak_channels[i]==0){
	    continue;
	  }
	  peak_number++;
	  found_baseline=matched_baseline=1;
	  REPORT("\n");
	  /* make a report about this peak */
	  strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(cycle_data->observation_time)));
	  baseline_to_antennae(cycle_data->baseline,&ant1,&ant2);
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "PEAK");
	  snprintf(tmp,BUFSIZE,"%04d",peak_number);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "number",BAD_CAST tmp);
	  REPORT("=P%04d= Observing source %s at %s, found peak in IF %d on longest baseline (ants %d & %d), in lag channel %d\n",
		 peak_number,cycle_data->sourcename,timestring,cycle_data->if_no,ant1,ant2,
		 cycle_data->spectrum_info.peak_channels[i]);
	  REPORT("=P%04d=  Peak details: amplitude = %.2f; spectrum average amplitude = %.2f, RMS noise level = %.2f, significance = %.2f sigma\n",
		 peak_number,cycle_data->pols_combined[cycle_data->spectrum_info.peak_channels[i]],
		 cycle_data->spectrum_info.estimated_average_amplitude,
		 cycle_data->spectrum_info.estimated_rms_noise,
		 (cycle_data->pols_combined[cycle_data->spectrum_info.peak_channels[i]]-
		  cycle_data->spectrum_info.estimated_average_amplitude)/
		 cycle_data->spectrum_info.estimated_rms_noise);
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "FOUND");
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "source_name",BAD_CAST cycle_data->sourcename);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "observation_time",BAD_CAST timestring);
	  snprintf(tmp,BUFSIZE,"%d",cycle_data->if_no);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "if_number",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%d",ant1);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "baseline_ant1",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%d",ant2);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "baseline_ant2",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",sqrtf(powf(cycle_data->u,2)+
					    powf(cycle_data->v,2)+
					    powf(cycle_data->w,2)));
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "baseline_length",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%d",cycle_data->spectrum_info.peak_channels[i]);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "lag_channel",BAD_CAST tmp);
	  /* what is the centre lag for the channel on the longest baseline? */
	  lag_calculator(cycle_data->spectrum_info.peak_channels[i],
			 cycle_data->if_bandwidth,&long_baseline_lag,
			 &long_baseline_lag_err);
/* 	  long_baseline_lag=(2*cycle_data->spectrum_info.peak_channels[i]+1)/ */
/* 	    (2*cycle_data->if_bandwidth); */
/* 	  long_baseline_lag_err=1/(2*cycle_data->if_bandwidth); */
	  snprintf(tmp,BUFSIZE,"%.2f",long_baseline_lag*1E9);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "baseline_lag",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",long_baseline_lag_err*1E9);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "baseline_lag_err",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",cycle_data->pols_combined[cycle_data->spectrum_info.peak_channels[i]]);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "amplitude",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",cycle_data->spectrum_info.estimated_average_amplitude);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "spectrum_average_amplitude",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",cycle_data->spectrum_info.estimated_rms_noise);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "spectrum_rms_noise",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.2f",(cycle_data->pols_combined[cycle_data->spectrum_info.peak_channels[i]]-
				       cycle_data->spectrum_info.estimated_average_amplitude)/
		   cycle_data->spectrum_info.estimated_rms_noise);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
	      BAD_CAST "significance",BAD_CAST tmp);
	  /* now check whether the peak is present in both polarisations */
	  additional_peak_checks(cycle_data,i);
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	  cycle_data->spectrum_info.peak_number[i]=peak_number;
	  /* check that this peak appears in the other IF at the same lag channel
	     (since lags are frequency independent) */
	  found_corresponding_peak=0;
	  for (second_cycle=*all_transforms;second_cycle!=NULL;
	       second_cycle=second_cycle->next){
	    if ((second_cycle->baseline==cycle_data->baseline)&&
		(second_cycle->observation_time==cycle_data->observation_time)&&
		(second_cycle->if_no==2)){
	      if (second_cycle->spectrum_info.npeaks_found>0){
		/* well there is a peak here at least */
		for (j=0;j<second_cycle->spectrum_info.npeaks_found;j++){
		  if (second_cycle->spectrum_info.peak_channels[j]==
		      cycle_data->spectrum_info.peak_channels[i]){
		    /* there is a peak in the same lag channel, that's good */
		    found_baseline++;
		    matched_baseline++;
		    strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(second_cycle->observation_time)));
		    baseline_to_antennae(second_cycle->baseline,&ant1,&ant2);
		    REPORT("=P%04d= Observing source %s at %s, found peak in IF %d on longest baseline (ants %d & %d), in lag channel %d\n",
			   peak_number,second_cycle->sourcename,timestring,second_cycle->if_no,ant1,ant2,
			   second_cycle->spectrum_info.peak_channels[j]);
		    REPORT("=P%04d=  Peak details: amplitude = %.2f; spectrum average amplitude = %.2f, RMS noise level = %.2f, significance = %.2f sigma\n",
			   peak_number,
			   second_cycle->pols_combined[second_cycle->spectrum_info.peak_channels[j]],
			   second_cycle->spectrum_info.estimated_average_amplitude,
			   second_cycle->spectrum_info.estimated_rms_noise,
			   (second_cycle->pols_combined[second_cycle->spectrum_info.peak_channels[j]]-
			    second_cycle->spectrum_info.estimated_average_amplitude)/
			   second_cycle->spectrum_info.estimated_rms_noise);
		    XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "FOUND");
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "source_name",BAD_CAST second_cycle->sourcename);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "observation_time",BAD_CAST timestring);
		    snprintf(tmp,BUFSIZE,"%d",second_cycle->if_no);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "if_number",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%d",ant1);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "baseline_ant1",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%d",ant2);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "baseline_ant2",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",sqrtf(powf(second_cycle->u,2)+
						      powf(second_cycle->v,2)+
						      powf(second_cycle->w,2)));
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "baseline_length",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%d",second_cycle->spectrum_info.peak_channels[j]);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "lag_channel",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",long_baseline_lag*1E9);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "baseline_lag",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",long_baseline_lag_err*1E9);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "baseline_lag_err",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",second_cycle->pols_combined[second_cycle->spectrum_info.peak_channels[j]]);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "amplitude",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",second_cycle->spectrum_info.estimated_average_amplitude);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "spectrum_average_amplitude",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",second_cycle->spectrum_info.estimated_rms_noise);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "spectrum_rms_noise",BAD_CAST tmp);
		    snprintf(tmp,BUFSIZE,"%.2f",(second_cycle->pols_combined[second_cycle->spectrum_info.peak_channels[j]]-
						 second_cycle->spectrum_info.estimated_average_amplitude)/
			     second_cycle->spectrum_info.estimated_rms_noise);
		    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
			      BAD_CAST "significance",BAD_CAST tmp);
		    /* now check whether the peak is present in both polarisations */
		    additional_peak_checks(cycle_data,i);
		    XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
		    second_cycle->spectrum_info.peak_number[j]=peak_number;
		    /* so now we need to check that the lag appears in the other baselines
		       at the correct lag */
		    second_cycle->spectrum_info.peak_verified[j]=PEAK_VERIFIED_FIRST_STAGE;
		    found_corresponding_peak=1;
		    long_baseline_length=sqrtf(powf(cycle_data->u,2)+
					       powf(cycle_data->v,2)+
					       powf(cycle_data->w,2));
		    for (baseline_cycle=*all_transforms;baseline_cycle!=NULL;
			 baseline_cycle=baseline_cycle->next){
		      if ((baseline_cycle->observation_time==cycle_data->observation_time)&&
			  (baseline_cycle->islongest_baseline!=1)){
			found_baseline++;
			/* figure out what the baseline ratio is */
			this_baseline_length=sqrtf(powf(baseline_cycle->u,2)+
						   powf(baseline_cycle->v,2)+
						   powf(baseline_cycle->w,2));
			baseline_ratio=this_baseline_length/long_baseline_length;
			REPORT("=P%04d= Scaling from long baseline length %.2f m to %.2f m on this baseline (ratio %.2f)\n",
			       peak_number,long_baseline_length,this_baseline_length,
			       baseline_ratio);
			/* now scale the lag found on the longest baseline and figure out
			   which lag channel the peak should appear in */
			lag_scaler(cycle_data->if_bandwidth,cycle_data->spectrum_info.peak_channels[i],
				   cycle_data->u,cycle_data->v,cycle_data->w,
				   baseline_cycle->u,baseline_cycle->v,baseline_cycle->w,
				   &this_baseline_expected_lag,&this_baseline_expected_lag_err,NULL);
/* 			this_baseline_expected_lag=this_baseline_length/long_baseline_length* */
/* 			  long_baseline_lag; */
/* 			this_baseline_expected_lag_err=this_baseline_expected_lag*long_baseline_lag_err/ */
/* 			  long_baseline_lag; */
			REPORT("=P%04d=  Long baseline lag is %.2f +/- %.2f ns, this baseline lag should be %.2f +/- %.2f ns\n",
			       peak_number,long_baseline_lag*1E9,long_baseline_lag_err*1E9,
			       this_baseline_expected_lag*1E9,this_baseline_expected_lag_err*1E9);
			/* now see if there are peaks on this baseline */
			found_peak_match=-1;
			if (baseline_cycle->spectrum_info.npeaks_found>0){
			  /* there is a peak, good start, now figure out what lag they are at */
			  for (k=0;k<baseline_cycle->spectrum_info.npeaks_found;k++){
			    lag_calculator(baseline_cycle->spectrum_info.peak_channels[k],
					   cycle_data->if_bandwidth,&this_baseline_lag,
					   &this_baseline_lag_err);
/* 			    this_baseline_lag=(2*baseline_cycle->spectrum_info.peak_channels[k]+1)/ */
/* 			      (2*cycle_data->if_bandwidth); */
/* 			    this_baseline_lag_err=1/(2*cycle_data->if_bandwidth); */
			    if ((check_for_overlap(this_baseline_lag,this_baseline_lag_err,
						   this_baseline_expected_lag,
						   this_baseline_expected_lag_err))==1){
			      /* we've found a matching peak! */
			      matched_baseline++;
			      strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",gmtime(&(baseline_cycle->observation_time)));
			      baseline_to_antennae(baseline_cycle->baseline,&ant1,&ant2);
			      REPORT("=P%04d=   Observing source %s at %s, found peak in IF %d on %.2f m baseline (ants %d & %d), in lag channel %d (%.2f +/- %.2f ns)\n",
				     peak_number,baseline_cycle->sourcename,timestring,baseline_cycle->if_no,
				     this_baseline_length,ant1,ant2,
				     baseline_cycle->spectrum_info.peak_channels[k],this_baseline_lag*1E9,
				     this_baseline_lag_err*1E9);
			      REPORT("=P%04d=    Peak details: amplitude = %.2f; spectrum average amplitude = %.2f, RMS noise level = %.2f, significance = %.2f sigma\n",
				     peak_number,
				     baseline_cycle->pols_combined[baseline_cycle->spectrum_info.peak_channels[k]],
				     baseline_cycle->spectrum_info.estimated_average_amplitude,
				     baseline_cycle->spectrum_info.estimated_rms_noise,
				     (baseline_cycle->pols_combined[baseline_cycle->spectrum_info.peak_channels[k]]-
				      baseline_cycle->spectrum_info.estimated_average_amplitude)/
				     baseline_cycle->spectrum_info.estimated_rms_noise);
			      XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "FOUND");
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "source_name",BAD_CAST baseline_cycle->sourcename);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "observation_time",BAD_CAST timestring);
			      snprintf(tmp,BUFSIZE,"%d",baseline_cycle->if_no);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "if_number",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%d",ant1);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_ant1",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%d",ant2);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_ant2",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",sqrtf(powf(baseline_cycle->u,2)+
								powf(baseline_cycle->v,2)+
								powf(baseline_cycle->w,2)));
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_length",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",this_baseline_expected_lag*1E9);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_expected_lag",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",this_baseline_expected_lag_err*1E9);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_expected_lag_err",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%d",baseline_cycle->spectrum_info.peak_channels[k]);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "lag_channel",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",this_baseline_lag*1E9);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_lag",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",this_baseline_lag_err*1E9);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "baseline_lag_err",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",baseline_cycle->pols_combined[baseline_cycle->spectrum_info.peak_channels[k]]);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "amplitude",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",baseline_cycle->spectrum_info.estimated_average_amplitude);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "spectrum_average_amplitude",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",baseline_cycle->spectrum_info.estimated_rms_noise);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "spectrum_rms_noise",BAD_CAST tmp);
			      snprintf(tmp,BUFSIZE,"%.2f",(baseline_cycle->pols_combined[baseline_cycle->spectrum_info.peak_channels[k]]-
							   baseline_cycle->spectrum_info.estimated_average_amplitude)/
				       baseline_cycle->spectrum_info.estimated_rms_noise);
			      XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
					BAD_CAST "significance",BAD_CAST tmp);
			      /* now check whether the peak is present in both polarisations */
			      additional_peak_checks(cycle_data,i);
			      XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
			      found_peak_match=k;
			      baseline_cycle->spectrum_info.peak_number[k]=peak_number;
			      break;
			    }
			  }
			}
			if (found_peak_match>-1){
			  /* we have found a peak, so mark the channel as a possibility */
			  baseline_cycle->spectrum_info.peak_verified[found_peak_match]=
			    PEAK_VERIFIED_FIRST_STAGE;
			}
		      }
		    }
		  }
		}
	      }
	      break;
	    }
	  }
	  if ((((float)matched_baseline/(float)found_baseline)>=arguments.minimum_matching_fraction)&&
	      (found_baseline>2))
	    /* this is good enough to call this peak a possible transient */
	    found_corresponding_peak=1;
	  else
	    found_corresponding_peak=0;
	  REPORT("=P%04d= Found %d peaks on %d baselines (%.2f %%); ",
		 peak_number,matched_baseline,found_baseline,((float)matched_baseline/
							      (float)found_baseline)*100);
	  XMLREPORT(xmlTextWriterStartElement,arguments.xml_report_writer,BAD_CAST "SUMMARY");
	  snprintf(tmp,BUFSIZE,"%d",found_baseline);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "available_baselines",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%d",matched_baseline);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "peak_baselines",BAD_CAST tmp);
	  snprintf(tmp,BUFSIZE,"%.4f",(float)matched_baseline/(float)found_baseline);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST "peak_baselines_fraction",BAD_CAST tmp);
	  if (found_corresponding_peak==0){
	    REPORT("this is not a transient candidate\n");
	    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		      BAD_CAST "transient_candidate",BAD_CAST "no");
	    /* oh well, no peak, we'll have to remove the peak found in the
	       cycle_data spectrum */
	    for (j=i+1;j<cycle_data->spectrum_info.npeaks_found;j++){
	      /* shift the peaks found down in the array by one */
	      cycle_data->spectrum_info.peak_channels[j-1]=
		cycle_data->spectrum_info.peak_channels[j];
	      cycle_data->spectrum_info.peak_verified[j-1]=
		cycle_data->spectrum_info.peak_verified[j];
	      cycle_data->spectrum_info.peak_number[j-1]=
		cycle_data->spectrum_info.peak_number[j];
	      cycle_data->spectrum_info.peak_present_both_pols[j-1]=
		cycle_data->spectrum_info.peak_present_both_pols[j];
	    }
	    cycle_data->spectrum_info.npeaks_found--;
	    /* we also have to reset all the other peaks we marked as first stage
	       possibilities to not possible */
	    for (second_cycle=*all_transforms;second_cycle!=NULL;
		 second_cycle=second_cycle->next){
	      for (j=0;j<second_cycle->spectrum_info.npeaks_found;j++){
		if (second_cycle->spectrum_info.peak_verified[j]==PEAK_VERIFIED_FIRST_STAGE)
		  second_cycle->spectrum_info.peak_verified[j]=PEAK_VERIFIED_NO;
	      }
	    }
	    if (cycle_data->spectrum_info.npeaks_found>0){
	      REALLOC(cycle_data->spectrum_info.peak_channels,
		      cycle_data->spectrum_info.npeaks_found*sizeof(int));
	      REALLOC(cycle_data->spectrum_info.peak_verified,
		      cycle_data->spectrum_info.npeaks_found*sizeof(int));
	      REALLOC(cycle_data->spectrum_info.peak_number,
		      cycle_data->spectrum_info.npeaks_found*sizeof(int));
	      REALLOC(cycle_data->spectrum_info.peak_present_both_pols,
		      cycle_data->spectrum_info.npeaks_found*sizeof(int));
	    } else {
	      FREE(cycle_data->spectrum_info.peak_channels);
	      FREE(cycle_data->spectrum_info.peak_verified);
	      FREE(cycle_data->spectrum_info.peak_number);
	      FREE(cycle_data->spectrum_info.peak_present_both_pols);
	    }
	    i--; /* need to keep checking peaks! */
	  } else {
	    /* we have a peak that works! */
	    /* we mark all the peaks that we've found as verified */
	    REPORT("this is a transient candidate\n");
	    XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		      BAD_CAST "transient_candidate",BAD_CAST "yes");
	    cycle_data->spectrum_info.peak_verified[i]=PEAK_VERIFIED_YES;
	    for (second_cycle=*all_transforms;second_cycle!=NULL;
		 second_cycle=second_cycle->next){
	      for (j=0;j<second_cycle->spectrum_info.npeaks_found;j++){
		if (second_cycle->spectrum_info.peak_verified[j]==PEAK_VERIFIED_FIRST_STAGE)
		  second_cycle->spectrum_info.peak_verified[j]=PEAK_VERIFIED_YES;
	      }
	    }
	  }
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	  XMLREPORT(xmlTextWriterEndElement,arguments.xml_report_writer);
	}
      }
    }
    cycle_data=cycle_data->next;
  }

  /* now go through the entire list again and remove all peaks that aren't 
     fully verified */
  for (cycle_data=*all_transforms;cycle_data!=NULL;cycle_data=cycle_data->next){
    for (j=0;j<cycle_data->spectrum_info.npeaks_found;j++){
      if (cycle_data->spectrum_info.peak_verified[j]!=PEAK_VERIFIED_YES){
	/* get rid of it */
	for (k=j+1;k<cycle_data->spectrum_info.npeaks_found;k++){
	  /* shift the peaks found down in the array by one */
	  cycle_data->spectrum_info.peak_channels[k-1]=
	    cycle_data->spectrum_info.peak_channels[k];
	  cycle_data->spectrum_info.peak_verified[k-1]=
	    cycle_data->spectrum_info.peak_verified[k];
	  cycle_data->spectrum_info.peak_number[k-1]=
	    cycle_data->spectrum_info.peak_number[k];
	  cycle_data->spectrum_info.peak_present_both_pols[k-1]=
	    cycle_data->spectrum_info.peak_present_both_pols[k];
	}
	cycle_data->spectrum_info.npeaks_found--;
	if (cycle_data->spectrum_info.npeaks_found>0){
	  REALLOC(cycle_data->spectrum_info.peak_channels,
		  cycle_data->spectrum_info.npeaks_found*sizeof(int));
	  REALLOC(cycle_data->spectrum_info.peak_verified,
		  cycle_data->spectrum_info.npeaks_found*sizeof(int));
	  REALLOC(cycle_data->spectrum_info.peak_number,
		  cycle_data->spectrum_info.npeaks_found*sizeof(int));
	  REALLOC(cycle_data->spectrum_info.peak_present_both_pols,
		  cycle_data->spectrum_info.npeaks_found*sizeof(int));
	} else {
	  FREE(cycle_data->spectrum_info.peak_channels);
	  FREE(cycle_data->spectrum_info.peak_verified);
	  FREE(cycle_data->spectrum_info.peak_number);
	  FREE(cycle_data->spectrum_info.peak_present_both_pols);
	}
	j--;
      }
    }
  }
  
  return(0);
}

int check_for_overlap(float value1,float value1_err,float value2,float value2_err){
  /* this routine figures out whether value1 and value2 are within their relative
     uncertainties of each other, returns 1 if they are, 0 otherwise */
  float min_value1,max_value1,min_value2,max_value2;
  int is_inrange=0;

  min_value1=value1-fabsf(value1_err);
  max_value1=value1+fabsf(value1_err);
  min_value2=value2-fabsf(value2_err);
  max_value2=value2+fabsf(value2_err);

  if ((min_value1>=min_value2)&&(min_value1<=max_value2))
    is_inrange=1;
  if ((value1>=min_value2)&&(value1<=max_value2))
    is_inrange=1;
  if ((max_value1>=min_value2)&&(max_value1<=max_value2))
    is_inrange=1;

  return(is_inrange);
}


int time_sort_data(correlation_data **data){
  /* this routine sorts the data structure and puts the earliest observation 
     at the top of the list, and the latest observation at the end of the list.
     it does this using the mergesort algorithm */
  correlation_data *p=NULL,*q=NULL,*e=NULL,*tail=NULL;
  int insize,nmerges,psize,qsize,i;

  if (*data==NULL)
    return(0);

  insize=1;

  while(1){
    p=*data;
    *data=NULL;
    tail=NULL;

    nmerges=0; /* count number of merges we do in this pass */
    
    while(p){
      nmerges++; /* there exists a merge to be done */
      /* step `insize' places along from p */
      q=p;
      psize=0;
      for (i=0;i<insize;i++){
	psize++;
	q=q->next;
	if (!q) break;
      }

      /* if q hasn't fallen off end, we have two lists to merge */
      qsize=insize;

      /* now we have two lists; merge them */
      while ((psize>0) || ((qsize>0) && q)){
	/* decide whether next element of merge comes from p or q */
	if (psize==0){
	  /* p is empty; e must come from q */
	  e=q;
	  q=q->next;
	  qsize--;
	} else if ((qsize==0) || !q){
	  /* q is empty; e must come from p */
	  e=p;
	  p=p->next;
	  psize--;
	} else if ((p->observation_time-q->observation_time)<=0){
	  /* first element of p is lower (or same); 
	     e must come from p */
	  e=p;
	  p=p->next;
	  psize--;
	} else {
	  /* first element of q is lower; e must come from q */
	  e=q;
	  q=q->next;
	  qsize--;
	}

	/* add the next element to the merged list */
	if (tail){
	  tail->next=e;
	} else {
	  *data=e;
	}
	/* maintain reverse pointers in a doubly linked list */
	e->prev=tail;
	tail=e;
      }

      /* now p has stepped `insize' places along, and q has too */
      p=q;
    }
    tail->next=NULL;

    /* if we have done only one merge, we're finished */
    if (nmerges<=1) /* all for nmerges==0, the empty list case */
      return(0);

    /* otherwise repeat, merging lists twice the size */
    insize*=2;
  }

  return(0);
}

float standard_deviation(float **values,int nelements,float *average){
  /* this routine calculates the standard deviation of the array passed to it */
  float value_stddev=0.0;
  int i;

  /* first get the mean value of the array */
  *average=meanvalue(values,nelements);

  /* now take the sum of each point minus the average */
  for (i=1;i<nelements;i++)
    /* DO NOT USE THIS ROUTINE AS A GENERAL PURPOSE STANDARD DEVIATION FINDER
       it ignores channel 0, since it will often be contaminated with a large peak
       whenever there is a continuum source at the phase centre of the observation */
    value_stddev+=powf(((*values)[i]-*average),2);
  /* divide by the number of elements */
  value_stddev/=nelements;
  /* and take the square root */
  value_stddev=sqrtf(value_stddev);

  return(value_stddev);
}

float meanvalue(float **values,int nelements){
  /* this routine calculates the average value of the array passed to it */
  float value_mean=0.0;
  int i;
  for (i=1;i<nelements;i++)
    /* DO NOT USE THIS ROUTINE AS A GENERAL PURPOSE AVERAGE FINDER
       it ignores channel 0, since it will often be contaminated with a large peak
       whenever there is a continuum source at the phase centre of the observation */
    value_mean+=(*values)[i];
  value_mean/=nelements;
  return(value_mean);
}

float standard_deviation_exclude(float **values,int nelements,float *average,
				 int n_exclusions,int **exclusions){
  /* this routine takes out some channels from the values array and then
     gets the standard deviation and mean of the data, so we can exclude
     known peak channels from the calculation */
  int i,j,k,new_size,y;
  float *new_values=NULL,retval;

  /* we need to check that channel 0 is not listed as an excluded
     channel, as we will have to include that channel ourselves */
  new_size=nelements-n_exclusions;
  for (i=0;i<n_exclusions;i++){
    if ((*exclusions)[i]==0)
      new_size++;
  }

  MALLOC(new_values,new_size*sizeof(float));
  for (i=0,j=0;i<nelements;i++){
    y=1;
    for (k=0;k<n_exclusions;k++){
      if (i==(*exclusions)[k])
	y=0;
    }
    if ((y==1)||(i==0)){
      /* note how we always include channel 0, as the standard deviation
         routine will automatically exclude that channel, so we need to
         give it some rubbish value anyway */
      new_values[j]=(*values)[i];
      j++;
    }
  }
  
  /* now send it to the normal routines */
  retval=standard_deviation(&new_values,new_size,average);
  FREE(new_values);
  
  /* and return */
  return(retval);
}

archive_file *file_info(char *filename){
  archive_file *cycle_file=NULL;

  for (cycle_file=arguments.existing_xml_report_files;cycle_file!=NULL;
       cycle_file=cycle_file->next){
    if (strstr(filename,cycle_file->filename)!=NULL)
      return(cycle_file);
  }
  return(NULL);
}

void build_file_database(void){
  xmlTextReaderPtr reader;
  xmlChar *name=NULL,*attrname=NULL,*attrvalue=NULL;
  int ret,nodetype,slash_offset,depth;
  archive_file *new_file=NULL;
  transient_peak *new_peak=NULL;
  struct tm observation_time;
  peak_baseline *new_baseline=NULL;

  reader=xmlReaderForFile(arguments.existing_xml_report,"ISO-8859-1",
			  XML_PARSE_RECOVER|XML_PARSE_NOENT|XML_PARSE_XINCLUDE);
  if (reader==NULL){
    ERROR("%s: cannot open file\n",arguments.existing_xml_report);
    arguments.use_existing_xml_report=0;
    return;
  }

  while((ret=xmlTextReaderRead(reader))==1){
    name=xmlTextReaderName(reader);
    nodetype=xmlTextReaderNodeType(reader);
    depth=xmlTextReaderDepth(reader);
    if (strcmp((char*)name,"FILE")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	/* the start of a new file */
	/* make a new file structure */
	MALLOC(new_file,sizeof(archive_file));
	new_file->transient_peaks=NULL;
	new_file->next=NULL;
	new_file->n_peaks=0;
	new_file->plot_name[0]='\0';
	/* get the name of the file from the attributes */
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  if (strcmp((char*)attrname,"filename")==0){
	    attrvalue=xmlTextReaderValue(reader);
	    if (attrvalue!=NULL){
	      for (slash_offset=strlen((char*)attrvalue);
		   slash_offset>0;slash_offset--){
		if (attrvalue[slash_offset]=='/'){
		  break;
		}
	      }
	      if (new_file!=NULL){
		strncpy(new_file->path,(char*)attrvalue,slash_offset);
		strcpy(new_file->filename,(char*)attrvalue+slash_offset+1);
	      }
	      xmlFree(attrvalue);
	    }
	  }
	  xmlFree(attrname);
	}
      } else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	/* the end of a file */
	/* add the file to the head of the list */
	if (new_file!=NULL){
	  new_file->next=arguments.existing_xml_report_files;
	  arguments.existing_xml_report_files=new_file;
	}
	new_file=NULL;
      }
    } else if (strcmp((char*)name,"SUMMARY")==0){
      if (depth==2){
	/* a summary of if there is any usable data in the file */
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  if (strcmp((char*)attrname,"usable")==0){
	    attrvalue=xmlTextReaderValue(reader);
	    if (attrvalue!=NULL){
	      if (strcmp((char*)attrvalue,"yes")==0){
		if (new_file!=NULL)
		  new_file->usable=YES;
	      } else if (strcmp((char*)attrvalue,"no")==0){
		if (new_file!=NULL)
		  new_file->usable=NO;
	      }
	      xmlFree(attrvalue);
	    }
	  }
	  xmlFree(attrname);
	}
      } else if (depth==3){
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  if (new_peak!=NULL){
	    attrname=xmlTextReaderName(reader);
	    attrvalue=xmlTextReaderValue(reader);
	    if (strcmp((char*)attrname,"available_baselines")==0){
	      new_peak->available_baselines=atoi((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"peak_baselines")==0){
	      new_peak->peak_baselines=atoi((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"peak_baselines_fraction")==0){
	      new_peak->peak_baselines_fraction=atof((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"transient_candidate")==0){
	      if (strcmp((char*)attrvalue,"yes")==0){
		new_peak->transient_candidate=YES;
	      } else if (strcmp((char*)attrvalue,"no")==0){
		new_peak->transient_candidate=NO;
	      }
	    }
	  }
	}
      }
    } else if (strcmp((char*)name,"PEAK")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	/* a peak number */
	MALLOC(new_peak,sizeof(transient_peak));
	new_peak->n_lags=0;
	new_peak->lag_peaks=NULL;
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (attrvalue==NULL){
	    /* oops this is a bad thing! */
	    xmlFree(attrname);
	    free(new_peak);
	    new_peak=NULL;
	    break;
	  }
	  if (strcmp((char*)attrname,"number")==0){
	    new_peak->peak_number=atoi((char*)attrvalue);
	  }
	}
      } else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	/* attach the peak to the file structure */
	if (new_peak!=NULL){
	  new_file->n_peaks++;
	  REALLOC(new_file->transient_peaks,
		  new_file->n_peaks*sizeof(transient_peak));
	  new_file->transient_peaks[new_file->n_peaks-1]=new_peak;
	  new_peak=NULL;
	}
      }
    } else if (strcmp((char*)name,"FOUND")==0){
      /* we've got a peak */
      if (new_peak!=NULL){
	MALLOC(new_baseline,sizeof(peak_baseline));
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (attrvalue==NULL){
	    /* oops this is a bad thing! */
	    xmlFree(attrname);
	    FREE(new_peak);
	    new_peak=NULL;
	    break;
	  }
	  if (strcmp((char*)attrname,"source_name")==0){
	    strcpy(new_peak->source_name,(char*)attrvalue);
	  } else if (strcmp((char*)attrname,"observation_time")==0){
	    sscanf((char*)attrvalue,"%2d:%2d:%2d %2d/%2d/%4d",
		   &(observation_time.tm_hour),&(observation_time.tm_min),
		   &(observation_time.tm_sec),&(observation_time.tm_mday),
		   &(observation_time.tm_mon),&(observation_time.tm_year));
	    observation_time.tm_mon--;
	    observation_time.tm_year-=1900;
	    new_peak->observation_time=timegm(&observation_time);
	  } else if (strcmp((char*)attrname,"if_number")==0){
	    new_baseline->if_number=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_ant1")==0){
	    new_baseline->baseline_ant1=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_ant2")==0){
	    new_baseline->baseline_ant2=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_length")==0){
	    new_baseline->baseline_length=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"lag_channel")==0){
	    new_baseline->lag_channel=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_lag")==0){
	    new_baseline->baseline_lag=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_lag_err")==0){
	    new_baseline->baseline_lag_err=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"amplitude")==0){
	    new_baseline->amplitude=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"spectrum_average_amplitude")==0){
	    new_baseline->spectrum_average_amplitude=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"spectrum_rms_noise")==0){
	    new_baseline->spectrum_rms_noise=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"significance")==0){
	    new_baseline->significance=atof((char*)attrvalue);
	    new_peak->n_lags++;
	    REALLOC(new_peak->lag_peaks,
		    new_peak->n_lags*sizeof(peak_baseline));
	    new_peak->lag_peaks[new_peak->n_lags-1]=new_baseline;
	    new_baseline=NULL;
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
      }
    }
  }
  xmlFreeTextReader(reader);
}

void free_archive_file (archive_file **unwanted_file) {
  int i,j;
  for (i=0;i<(*unwanted_file)->n_peaks;i++){
    for (j=0;j<(*unwanted_file)->transient_peaks[i]->n_lags;j++){
      FREE((*unwanted_file)->transient_peaks[i]->lag_peaks[j]);
    }
    FREE((*unwanted_file)->transient_peaks[i]->lag_peaks);
  }
  FREE((*unwanted_file)->transient_peaks);
  FREE(*unwanted_file);
}

/**
 * ConvertInput:
 * @in: string in a given encoding
 * @encoding: the encoding used
 *
 * Converts @in into UTF-8 for processing with libxml2 APIs
 *
 * Returns the converted UTF-8 string, or NULL in case of error.
 */
xmlChar *
ConvertInput(const char *in, const char *encoding)
{
    xmlChar *out;
    int ret;
    int size;
    int out_size;
    int temp;
    xmlCharEncodingHandlerPtr handler;

    if (in == 0)
        return 0;

    handler = xmlFindCharEncodingHandler(encoding);

    if (!handler) {
        printf("ConvertInput: no encoding handler found for '%s'\n",
               encoding ? encoding : "");
        return 0;
    }

    size = (int) strlen(in) + 1;
    out_size = size * 2 - 1;
    out = (unsigned char *) xmlMalloc((size_t) out_size);

    if (out != 0) {
        temp = size - 1;
        ret = handler->input(out, &out_size, (const xmlChar *) in, &temp);
        if ((ret < 0) || (temp - size + 1)) {
            if (ret < 0) {
                printf("ConvertInput: conversion wasn't successful.\n");
            } else {
                printf
                    ("ConvertInput: conversion wasn't successful. converted: %i octets.\n",
                     temp);
            }

            xmlFree(out);
            out = 0;
        } else {
            out = (unsigned char *) xmlRealloc(out, out_size + 1);
            out[out_size] = 0;  /*null terminating out */
        }
    } else {
        printf("ConvertInput: no mem\n");
    }

    return out;
}

void lag_calculator(int channel_number,float bandwidth,float *lag,float *lag_err){
  /* this routine calculates the lag of a specified channel where the sampling
     rate is determined by the bandwidth in the usual Nyquist way */
  *lag=(2*(float)channel_number+(float)1)/(2*bandwidth);
  *lag_err=1/(2*bandwidth);
}

int channel_calculator(float bandwidth,float lag_value,float lag_err){
  /* this routine takes a lag and its associated error, and figures out which
     channel that lag corresponds to given the bandwidth of the observation */
  int lag_channel=-1,current_channel;
  float central_lag,lag_width;

  for (current_channel=0;current_channel<16384;current_channel++){
    lag_calculator(current_channel,bandwidth,&central_lag,&lag_width);
    if ((check_for_overlap(lag_value,lag_err,central_lag,lag_width))==1){
      /* this channel would contain this lag */
      lag_channel=current_channel;
      break;
    }
  }

  return(lag_channel);
}

void lag_scaler(float bandwidth,int long_baseline_channel,
		float long_baseline_u,float long_baseline_v,float long_baseline_w,
		float this_baseline_u,float this_baseline_v,float this_baseline_w,
		float *this_baseline_lag,float *this_baseline_lag_err,
		int *this_baseline_channel){
  /* this routine figures out what channel a peak should be in on any
     particular baseline, given knowledge of where the peak is on another
     baseline */
  float long_baseline_length,this_baseline_length,baseline_ratio;
  float long_baseline_lag,long_baseline_lag_err;
  float tmp_baseline_lag,tmp_baseline_lag_err;

  /* figure out the baseline lengths and the ratio between them */
  long_baseline_length=sqrtf(powf(long_baseline_u,2)+
			     powf(long_baseline_v,2)+
			     powf(long_baseline_w,2));
  this_baseline_length=sqrtf(powf(this_baseline_u,2)+
			     powf(this_baseline_v,2)+
			     powf(this_baseline_w,2));
  baseline_ratio=this_baseline_length/long_baseline_length;

  /* get the lag on the long baseline */
  lag_calculator(long_baseline_channel,bandwidth,&long_baseline_lag,
		 &long_baseline_lag_err);

  /* now scale the lag found on the longest baseline and figure out
     which lag channel the peak should appear in */
  /* we check for NULL pointers before we assign these values in case
     the calling routine doesn't actually care about them */
  tmp_baseline_lag=this_baseline_length/long_baseline_length*
    long_baseline_lag;
  tmp_baseline_lag_err=tmp_baseline_lag*long_baseline_lag_err/
    long_baseline_lag;
  if (this_baseline_lag!=NULL)
    *this_baseline_lag=tmp_baseline_lag;
  if (this_baseline_lag_err!=NULL)
    *this_baseline_lag_err=tmp_baseline_lag_err;

  /* which channel does this correspond to? */
  if (this_baseline_channel!=NULL)
    *this_baseline_channel=channel_calculator(bandwidth,tmp_baseline_lag,
					      tmp_baseline_lag_err);

}
 
void additional_peak_checks(correlation_data *transform,int peak_number){
  /* this routine performs some additional checks on the remaining
     peaks, removing any that fail the checks */
  correlation_data *cycle_data=NULL;
  int i,*excluded_channels=NULL;
  float rms_noise,average_value;
  char tmp[BUFSIZE],tmp2[BUFSIZE];

  if (transform==NULL)
    return;

  cycle_data=transform;
  /* the first thing we will check is that the peaks are found in both
     of the individual polarisations */
  if (cycle_data->spectrum_info.npeaks_found>0){
    /* there are peaks in this spectrum */
    MALLOC(excluded_channels,
	   cycle_data->spectrum_info.npeaks_found*sizeof(int));
    for (i=0;i<cycle_data->spectrum_info.npeaks_found;i++){
      /* add the peaks to the channels to exclude from the rms
	 noise calculation */
      excluded_channels[i]=cycle_data->spectrum_info.peak_channels[i];
    }
    for (i=0;i<cycle_data->if_nstok;i++){
      /* get the average level and rms noise level from the spectrum */
      if (cycle_data->individual_pols[i]!=NULL){
	rms_noise=standard_deviation_exclude(&(cycle_data->individual_pols[i]),
					     cycle_data->if_nfreq,
					     &average_value,
					     cycle_data->spectrum_info.npeaks_found,
					     &excluded_channels);
	snprintf(tmp,BUFSIZE,"%.2f",average_value);
	snprintf(tmp2,BUFSIZE,"polarisation_%s_average_amplitude",
		 cycle_data->stokes_name[i]);
	XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		  BAD_CAST tmp2,BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",rms_noise);
	snprintf(tmp2,BUFSIZE,"polarisation_%s_rms_noise",
		 cycle_data->stokes_name[i]);
	XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		  BAD_CAST tmp2,BAD_CAST tmp);
	/* now test to see whether the peak channels are well above the
	   average level in each polarisation (we will halve the significance
	   requirement here) */
	if (cycle_data->individual_pols[i][cycle_data->spectrum_info.peak_channels[peak_number]]<
	    (average_value+(arguments.minimum_peak_sigma_factor/2)*rms_noise)){
	  /* we failed the test, there is no significant peak here */
	  cycle_data->spectrum_info.peak_present_both_pols[peak_number]=0;
	  snprintf(tmp2,BUFSIZE,"polarisation_%s_detected",
		   cycle_data->stokes_name[i]);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST tmp2,BAD_CAST "no");
	} else {
	  snprintf(tmp2,BUFSIZE,"polarisation_%s_detected",
		   cycle_data->stokes_name[i]);
	  XMLREPORT(xmlTextWriterWriteAttribute,arguments.xml_report_writer,
		    BAD_CAST tmp2,BAD_CAST "yes");
	}
      }
    }
  }
}
